1function void boot_printout(int type, int n, int rep, 2 int failed, const matrix Spar_mat) 3 # Sep 2020: change from pointerized Spar_mat to const 4 5 matrix bm = meanc(Spar_mat) 6 matrix bV = mcov(Spar_mat) 7 scalar nc = cols(Spar_mat) 8 scalar n2 = n*n 9 10 # force numerical zeros to real zeros 11 e = diag(bV) .< 1.0e-12 12 if maxc(e) 13 e = selifc(seq(1, nc), e') 14 bV[e,] = 0 15 bV[,e] = 0 16 endif 17 18 printf "Bootstrap results (%d replications, %d failed)\n", \ 19 rep + failed, failed 20 21 if (type != 3) && ( nc == 2 * n2 ) # was cols(Spar_mat) 22 # Long-run matrix at the end exists! 23 # And so this is a model with long-run restrictions 24 # (Bl-Quah-style or SVEC), 25 # or the user forced it via calc_lr. 26 27 matrix bK = mshape( bm[1:n2], n, n) 28 printStrMat(bK, bV[1:n2, 1:n2], "C") 29 matrix bL = mshape( bm[n2+1:], n, n) 30 printStrMat(bL, bV[1+n2: , 1+n2: ], "LongRun") 31 32 elif type != 3 33 34 # C model without printing the long-run matrix 35 # (SVEC / type 4 should not happen here, because it has 36 # long-run constraints by construction) 37 38 matrix bK = mshape(bm,n,n) 39 printStrMat(bK, bV, "C") 40 41 elif type == 3 # AB model 42 matrix bmA = mshape( bm[1:n2], n, n) 43 printStrMat(bmA, bV[1:n2, 1:n2], "A") 44 matrix bmB = mshape( bm[n2+1:], n, n) 45 printStrMat(bmB, bV[1+n2:,1+n2:], "B") 46 47 else 48 funcerr "shouldn't happen" 49 endif 50end function 51 52############### 53 54function matrix bias_correction(const bundle b, 55 const matrix Y0, matrix *BC, 56 const bundle bparams[null]) 57 58 /* This function implements a bias correction for 59 the estimate of the VAR parameters as per Kilian, 60 REStat (1998). 61 (The SVEC case is not allowed here, but it must 62 be checked on the outside.) 63 64 Sep 2019: re-did interface to use bundle; the new 65 'boottype' member means: 66 1: standard resampling (default, as before) 67 2-4: wild bootstrap (3 variants) 68 5: moving-blocks bootstrap 69 70 (So the check below for what we need in the bundle.) 71 */ 72 73 errorif( !inbundle(b,"BCiter") || !inbundle(b,"VARpar") || \ 74 !inbundle(b,"E") || !inbundle(b,"X") || !inbundle(b,"boottype") || \ 75 !inbundle(b,"bmu"), "needed input missing in bundle arg") 76 77 # check for stationarity first 78 scalar maxmod = max_eval(b.VARpar) 79 80 if maxmod < 0.9999 81 82 matrix innov = prepres(b.U, b.boottype) 83 bundle bootstuff = _(Y0, innov) # was: defbundle("Y0",Y0, "innov",innov) 84 if exists(bparams) 85 bootstuff.bparams = bparams 86 endif 87 matrix Ab = avg_VARpar_boot(b, bootstuff)[1] 88 matrix BC = b.VARpar - Ab 89 add_and_smash(&Ab, BC) # was H = ..., unused 90 91 else # not stationary 92 matrix Ab = b.VARpar 93 endif 94 return Ab 95end function 96 97#################### 98 99function matrices avg_VARpar_boot(const bundle b, const bundle stuff) 100 # In principle we could call base_est() below in the loop, 101 # because that's what it's for, but since we only need 102 # the reduced-form AR coefficients, it's not worth it... 103 104 # This function returns an array (matrices) just in case in the 105 # future we want to apply the parallel_func MPI approach to it. 106 # (This is a requirement there.) 107 # Of course in essence it's still just a single matrix. 108 109 bundle bparams = inbundle(stuff, "bparams") ? stuff.bparams : null 110 111 # simulate the average reduced-form coeffs 112 matrix Absum = zeros(b.n, b.n*b.p) 113 matrix lagseq = seq(1,b.p) 114 matrix Uinit = zeros(b.p, b.n) 115 116 if rows(stuff.Y0) != b.p 117 funcerr "wrong assumption about length of Y initvals here!" 118 endif 119 120 loop i = 1..b.BCiter 121 matrix U = Uinit | drawbootres(stuff.innov, b.boottype, bparams) 122 if cols(b.bmu) > 0 # was mu 123 U += b.bmu 124 endif 125 matrix bY = varsimul(b.VARpar, U[b.p + 1: ,], stuff.Y0) # was rows(stuff.Y0)+1 126 matrix reg = b.X ~ mlag(bY, lagseq) 127 matrix Pi = mols(bY[b.p + 1: ,], reg[b.p + 1: ,]) 128 Absum += transp(Pi[b.k + 1: b.k + b.n*b.p, ]) 129 endloop 130 131 return defarray(Absum ./ b.BCiter) 132end function 133 134################## 135 136function void calc_bmu(bundle *obj) 137 # disentangle determ/exog 138 # Sep 2020: add the result directly to the (pointerized) bundle 139 140 matrix bmu = zeros(obj.T, obj.n) 141 if obj.k && obj.type == 4 # SVEC, 142 # (with some unrestr. exo apart from const/trend) 143 if obj.jcase == 1 144 bmu = obj.X * obj.mu 145 146 elif obj.jcase == 2 || obj.jcase == 3 147 bmu = obj.X * obj.mu[2:, ] 148 # need to add restr. or unrestr. const later 149 150 elif obj.jcase == 4 || obj.jcase == 5 151 bmu = obj.X * obj.mu[3:, ] 152 # need to add restr./unr. const & trend later 153 endif 154 155 elif obj.k # no SVEC 156 bmu = obj.X * obj.mu # this was the pre-1.4 handled case 157 endif 158 159 # more special treatment of SVEC 160 if obj.type == 4 161 # add constant 162 if obj.jcase > 2 # unrestricted 163 bmu = bmu .+ obj.mu[1, ] # (use broadcasting) 164 165 elif obj.jcase == 2 # restricted 166 bmu = bmu .+ (obj.jbeta[obj.n + 1, ] * obj.jalpha') 167 endif 168 169 # add trend 170 if obj.jcase == 4 # restricted 171 bmu += seq(1, obj.T)' obj.jbeta[obj.n + 1, ] * obj.jalpha' 172 173 elif obj.jcase == 5 # unrestricted 174 bmu += seq(1, obj.T)' obj.mu[2, ] 175 endif 176 endif 177 178 matrix obj.bmu = bmu 179end function 180 181#---------------------------------- 182# These private functions introduced in Aug/Sep 2019 183# for general bootstrap 184 185function matrix prepres(const matrix U, 186 int boottype[1:5:1] "bootstrap type" \ 187 {"resampling", "wildN", "wildR", "wildM", \ 188 "moving blocks"}) 189 # The input is demeaned (column-wise) only for 190 # the resampling (type 1). It isn't done together with the 191 # new draws to avoid repeating it many times. 192 193 # (Can now have 'const matrix U' as gretl versions >=2019d handle this fine 194 # when doing 'return U'.) 195 196 if boottype == 1 197 return cdemean(U) 198 else 199 return U 200 endif 201end function 202 203 204function matrix drawbootres(const matrix U, 205 int boottype[1:5:1] "bootstrap type" \ 206 {"resampling", "wildN", "wildR", "wildM", \ 207 "moving blocks"}, 208 const bundle bparams[null] "further inputs" ) 209 210 /* 211 Construct a new draw of residuals for bootstrapping, where U 212 is a Txn matrix. 213 U can be original residuals or can be some pre-processed 214 input. (See the prepres() function - the pre-processing is 215 not done here to avoid doing it repeatedly.) 216 217 Currently boottype can go from 1 to 5: 218 1: traditional residual resampling (using gretl's resample()) 219 2-4: wild bootstrap with the help of a standard normal, 220 Rademacher, Mammen distributions respectively 221 5: RBMBB Brüggemann/Jentsch/Trenkler 2016. A block length 222 must be chosen. 223 224 Other bootstrap types may be added in the future. 225 (The U input could be empty then, e.g. for a purely parametric- 226 distribution bootstrap. The bundle bparams is intended to be 227 used then to specify parameters, in a way that is to be determined.) 228 */ 229 230 # process extra input 231 bl = 0 # signals default data-based choice 232 if exists(bparams.moveblocklen) && boottype == 5 233 l = bparams.moveblocklen 234 endif 235 236 ## standard 237 if boottype == 1 238 return resample(U) 239 endif 240 241 T = rows(U) 242 243 ## wild 244 if boottype <= 4 245 246 if boottype == 2 247 # Normal 248 matrix w = mnormal(T) 249 elif boottype == 3 250 # Rademacher 251 matrix w = muniform(T) .< 0.5 ? 1 : -1 252 elif boottype == 4 253 # Mammen 254 scalar s5 = sqrt(5) 255 scalar p = (0.5/s5) * (s5 + 1) 256 matrix w = 0.5 + (muniform(T) .< p ? -s5/2 : s5/2) 257 endif 258 259 return U .* w 260 261 ## residual-based moving blocks 262 elif boottype == 5 263 if bl > T 264 print "Warning: block length too large, falling back to auto" 265 bl = 0 266 endif 267 if bl == 0 # use default 268 bl = xmax(2, floor(T/10)) 269 endif 270 271 # necessary number of blocks 272 s = ceil(T/bl) 273 # blocks starting points 274 matrix c = mrandgen(i, 1, T-bl+1, 1, s) 275 # convert to indices of all needed obs 276 matrix ndx = vec(c .+ seq(0, bl-1)') 277 278 ## recentring 279 matrix m = mshape(NA, bl, cols(U)) 280 # calculate respective averages (bl different ones) 281 loop i = 1..bl 282 m[i,] = meanc(U[i: T-bl+i,]) 283 endloop 284 285 ndx = ndx[1:T] # cut off "overhanging tail" 286 return U[ndx,] - (ones(s) ** m)[1:T, ] 287 288 else 289 funcerr "Shouldn't happen" 290 endif 291 292end function 293 294################## 295 296function void boottypechoice(bundle *mod, string btypestr) 297 temp = boottypecode(btypestr) 298 if temp == 0 299 print "Warning: ignoring unrecognized SVAR bootstrap type." 300 else 301 mod.boottype = temp 302 endif 303end function 304 305################# 306 307function void maybe_do_biascorr(bundle *bobj, const bundle bparams) 308 # this function adds "ABCorA" and perhaps "Psi" to bobj 309 310 if bobj.biascorr && (bobj.type != 4) # not available w/unit roots # BIASCORR 311 matrix Psi = {} 312 matrix start = bobj.Y[1: bobj.p, ] 313 matrix bobj.ABCorA = bias_correction(bobj, start, &Psi, bparams) # ABC # moved up from inside the loop 314 matrix bobj.Psi = Psi 315 else 316 matrix bobj.ABCorA = bobj.VARpar 317 endif 318 319end function 320 321/* ------------------------------------------------------------------- */ 322/* --- Main bootstrap function --------------------------------------- */ 323/* ------------------------------------------------------------------- */ 324 325function scalar SVAR_boot(bundle *obj, 326 int rep[0::2000] "bootstrap iterations", 327 scalar alpha[0:1:0.9] "CI coverage", 328 bool quiet[1], 329 string btypestr[null] "bootstrap type", 330 int biascorr[-1:2:-1] "bias correction (non-SVEC)") 331 332 # btypestr: can be "resample" / "resampling", 333 # "wildN"/"wild", "wildR", "wildM" 334 335 # The default value for biascorr of -1 means: 336 # Do not override the previous setting. 337 338 ## Copy some params and choices 339 loop foreach i n k T p type 340 scalar $i = obj.$i 341 endloop 342 obj.nboot = rep # record bootstrap details 343 obj.boot_alpha = alpha # into original model 344 345 errorif( type == 10, "Wrong turn: Set-ID not for bootstrapping...") 346 347 ## Bootstrap type choice (if different from default) 348 if exists(btypestr) 349 boottypechoice(&obj, btypestr) 350 endif 351 352 # Copy optional block length choice 353 bundle bparams = null 354 if obj.boottype == 5 && inbundle(obj, "moveblocklen") 355 bparams.moveblocklen = obj.moveblocklen 356 endif 357 358 # Bias correction choice, leave or update? 359 obj.biascorr = (biascorr == -1) ? obj.biascorr : biascorr 360 361 # define default for bias correction iterations 362 if obj.biascorr && !inbundle(obj, "BCiter") 363 obj.BCiter = 1024 364 endif 365 366 ## Various needed stuff 367 if type == 3 # AB 368 matrix bmA bmB # needed as memory for transfer 369 elif type == 4 370 matrix J = zeros(n - obj.crank, obj.crank) | I(obj.crank) 371 endif 372 matrix start = obj.Y[1:p, ] # Y0 373 374 # disentangle determ/exog: 375 calc_bmu(&obj) # adds obj.bmu 376 377 # store a copy of the model for bootstrap 378 bundle bobj = obj 379 380 # Do the bias correction pre-step if applicable 381 maybe_do_biascorr(&bobj, bparams) 382 383 printf "\nBootstrapping model (%d iterations)\n", rep 384 printf "Bootstrap type: %s\n", btypestring(obj.boottype) 385 printf "Bias correction: %s\n\n", BCstring(obj.biascorr) 386 flush 387 388 ## Actual bootstrap simulation 389 matrices bootout = SVAR_boot_innerloop(&bobj, obj, bparams) 390 /* 391 "bootirfs": 392 each bootstrap replication on one row; each row contains 393 the vectorisation of the complete IRF matrix 394 */ 395 matrix bootirfs = bootout[1] # zeros(rep, (h+1) * n2) 396 # Spar_mat is probably somewhat redundant..., only for the printout 397 matrix Spar_mat = bootout[2] # zeros(rep, n2) or zeros(rep, 2*n2) 398 scalar failed = bootout[3] 399 400 if !quiet 401 boot_printout(type, n, rep, failed, Spar_mat) 402 endif 403 404 # quantiles of bootstrapped IRFs used in graphs 405 q_alpha = 0.5 * (1 - alpha) # changed in v1.5 406 matrix locb = quantile(bootirfs, q_alpha) 407 matrix hicb = quantile(bootirfs, 1 - q_alpha) 408 matrix mdn = quantile(bootirfs, 0.5) 409 410 bundle bootdata = null 411 bootdata.rep = rep # no of replications 412 bootdata.alpha = alpha # alpha 413 bootdata.biascorr = obj.biascorr # type of bias correction 414 scalar h = obj.horizon 415 scalar n2 = n*n 416 matrix bootdata.lo_cb = mshape(locb, h+1, n2) # lower bounds 417 matrix bootdata.hi_cb = mshape(hicb, h+1, n2) # upper bounds 418 matrix bootdata.mdns = mshape(mdn, h+1, n2) # medians 419 420 bundle obj.bootdata = bootdata 421 return failed 422end function 423 424############################################ 425## SVAR boot inner loop ####### 426 427function matrices SVAR_boot_innerloop(bundle *bobj, const bundle obj, const bundle bparams) 428 429 # will return a three-matrix array: 430 # 1: bootirfs 431 # 2: Spar_mat 432 # 3: number of failures (1x1 scalar) 433 434 # copy 435 type = obj.type 436 n2 = obj.n * obj.n 437 438 # Prepare the residuals 439 matrix innov = prepres(obj.E, obj.boottype) 440 441 ## prepare output 442 matrix bootirfs = zeros(obj.nboot, (obj.horizon + 1) * n2) 443 # Spar_mat: the result matrix 444 # (-- type==10 (set-ID) cannot/should not happen in here --) 445 # need more cols if saving either A,B (for type 3) or C and the long-run matrix 446 numcols = type>2 || (rows(obj.Rd1l) || obj.calc_lr) ? 2*n2 : n2 447 matrix Spar_mat = zeros(obj.nboot, numcols) 448 449 if type == 3 && inbundle(bobj, "C") 450 # clean up object from pre-computed C matrix 451 delete bobj.C 452 endif 453 454 i = 1 455 failed = 0 456 set loop_maxiter 16384 457 loop while i <= obj.nboot 458 # clear previous bootstrap indicator 459 bobj.step = 0 460 461 # generate bootstrap disturbances (bmu may be zero) 462 matrix U = obj.bmu[obj.p + 1:, ] + drawbootres(innov, obj.boottype, bparams) 463 464 # generate bootstrap data and store it in bootstrap object 465 matrix bobj.Y = varsimul(bobj.ABCorA, U, obj.Y[1: obj.p, ] ) 466 467 # estimate VAR parameters, special treatment VECM/SVEC 468 if type == 4 469 vecm_est(&bobj) 470 else 471 base_est(&bobj) 472 endif 473 474 matrix bA = bobj.VARpar # estimates (first n rows of companion mat) 475 matrix bSigma = bobj.Sigma 476 matrix theta = obj.theta # init original SVAR params 477 # (C/A&B apparently, in suitable form...) 478 479 errcode = 0 480 481 ## Full bias correction 482 /* (The bc-ed VARpar need to be done before the new C matrix is 483 calculated at least if there are long-run constraints, because 484 then they enter via C1 through fullRd into C. 485 Otherwise the new C only depends on Sigma. - BTW, we do not update 486 Sigma in the full biascorr case. (In theory we could, by 487 re-calculating the residuals using the bc-ed VARpar. 488 We shouldn't, should we?) 489 */ 490 491 if obj.biascorr == 2 && type != 4 # only for non-SVEC 492 scalar H = add_and_smash(&bA, bobj.Psi) 493 494 if ok(H) 495 bobj.VARpar = bA 496 else 497 errcode = 101 498 endif 499 endif 500 501 /* now re-estimate C, according to model type */ 502 503 if type == 1 # Cholesky 504 matrix K = cholesky(bSigma) 505 506 elif type == 2 || type == 4 # consolidate the SVEC case here 507 508 # Watch out for the possibility of constraints that were 509 # partly redundant originally. 510 # 511 # In the AB-model the cleaned restrictions are inherited 512 # from the original model, apart from the fact that they 513 # should already be caught at specification time (in 514 # SVAR_restrict). 515 # 516 # But with long-run constraints (C, and especially SVEC 517 # with weakly exog.) the restrictions depend on the 518 # bootstrap data, so the cleaning of the restrictions must 519 # be re-done here. 520 521 if type != 3 && inbundle(obj, "cleanfullRd") 522 bobj.fullRd = {} # reset 523 id = ident(&bobj, 0) # re-creates fullRd and possibly cleanfullRd 524 if !id 525 printf "Ident problem in bootstrap draw %d\n", i 526 funcerr "Unexpected ID problem in bootstrap" 527 endif 528 529 elif type == 4 || nelem(obj.Rd1l) # long-run constraints 530 bobj.fullRd = get_full_Rd(&bobj, 0) # update fullRd 531 532 endif 533 534 if inbundle(bobj, "cleanfullRd") 535 matrix fullRd = bobj.cleanfullRd 536 else 537 matrix fullRd = bobj.fullRd 538 endif 539 540 matrix K = estC(&theta, bSigma, fullRd, null, &errcode, obj.optmeth, 0) 541 542 elif type == 3 # "AB" 543 matrices transferAB = array(0) 544 # matrix Rd2 = type==3 ? obj.Rd0 : obj.Rd1l 545 546 # new: get A,B instead of re-calc'ing it below 547 # (obj.Rd0 was Rd2 before, but redundant...) 548 matrix K = estAB(&theta, bSigma, obj.Rd0, obj.Rd1, null, \ 549 &errcode, obj.optmeth, 0, &transferAB) 550 endif 551 552 ## Process and store the simulated C results 553 if !errcode && rows(K) == obj.n 554 bobj.step = 2 555 bobj.theta = theta 556 557 # we don't treat the AB-model specially here (no reason to) 558 maybe_flip_columns(obj.C, &K) 559 560 if (type == 1) || (type == 2) || (type == 4) 561 bobj.C = K # is used in doIRF() 562 Spar_mat[i, 1: n2] = vec(K)' 563 564 /* New Oct 2017: Also bootstrap the long-run matrix if wanted 565 Jan 2018: add type 4 */ 566 if ( type < 3 && ( rows(bobj.Rd1l) || bobj.calc_lr ) ) || type == 4 567 # (a plain or C model w/ long-run constr, or user switch) 568 # long-run mat (C1 comes from get_full_Rd() above 569 # (except type 1)): 570 571 matrix C1 = (type == 2 || type == 4) ? \ 572 bobj.C1 : C1mat(bobj.VARpar) 573 matrix bobj.lrmat = C1 * bobj.C 574 575 # attach it to the other bootstrap result 576 Spar_mat[i, n2+1 : 2*n2] = vec(bobj.lrmat)' 577 endif 578 579 elif type == 3 580 # (Sven): the following stuff comes from estAB above 581 bobj.S1 = transferAB[1] 582 bobj.S2 = transferAB[2] 583 Spar_mat[i,] = vec(bobj.S1)' ~ vec(bobj.S2)' 584 endif 585 586 endif 587 588 if !errcode && rows(K) == obj.n 589 doIRF(&bobj) 590 bootirfs[i,] = vec(bobj.IRFs)' 591 i++ 592 else 593 failed++ 594 outfile stderr 595 printf "Iter %4d failed (error code = %d)\n", i, errcode 596 end outfile 597 endif 598 endloop 599 600 return defarray(bootirfs, Spar_mat, {failed}) 601end function 602