*! _mvardlurt_bootstrap — Bootstrap engine for mvardlurt *! Version 2.0.0 — 2026-02-24 *! *! Implements the parametric bootstrap from Sam, McNown, Goh & Goh (2024). *! Uses st_store/st_data approach — no tempfiles needed. *! Recursive DGP in Mata exactly replicates EViews model a / a.solve. capture program drop _mvardlurt_bootstrap program define _mvardlurt_bootstrap, rclass version 14 syntax varname(ts), /// INDEPvar(varname ts) /// independent variable Plag(integer) /// optimal dy lag length Qlag(integer) /// optimal dx lag length Case(integer) /// 1, 3, or 5 REPS(integer) /// bootstrap replications SEED(integer) /// random seed local depvar "`varlist'" local indepvar "`indepvar'" local opt_p = `plag' local opt_q = `qlag' // ========================================================================= // 1. SETUP // ========================================================================= qui tsset qui count local TT = r(N) set seed `seed' local nocons_opt "" if `case' == 1 local nocons_opt "noconstant" local opt_lag = `opt_p' + `opt_q' // Build regressor strings local dy_lags "" if `opt_p' > 0 { forvalues j = 1/`opt_p' { local dy_lags "`dy_lags' L`j'.D.`depvar'" } } local dx_lags "" if `opt_q' > 0 { forvalues j = 1/`opt_q' { local dx_lags "`dx_lags' L`j'.D.`indepvar'" } } local det_regs "" if `case' == 5 { tempvar ttrend qui gen double `ttrend' = _n local det_regs "`ttrend'" } // Model specifications (matching EViews exactly) local rest_t_reg "L.`indepvar' `dy_lags' `dx_lags' `det_regs'" local rest_f_reg "L.`depvar' `dy_lags' `dx_lags' `det_regs'" local full_reg "L.`depvar' L.`indepvar' `dy_lags' `dx_lags' `det_regs'" // ========================================================================= // 2. ESTIMATE RESTRICTED MODELS & EXTRACT/RECENTER RESIDUALS // ========================================================================= // t-restricted: impose y(-1) = 0 qui regress D.`depvar' `rest_t_reg', `nocons_opt' tempvar resid_t qui predict double `resid_t', residuals tempname b_rest_t mat `b_rest_t' = e(b) qui su `resid_t', meanonly qui replace `resid_t' = `resid_t' - r(mean) // f-restricted: impose x(-1) = 0 qui regress D.`depvar' `rest_f_reg', `nocons_opt' tempvar resid_f qui predict double `resid_f', residuals tempname b_rest_f mat `b_rest_f' = e(b) qui su `resid_f', meanonly qui replace `resid_f' = `resid_f' - r(mean) // ========================================================================= // 3. COPY EVERYTHING TO MATA (once, before the loop) // ========================================================================= // Find the column index of depvar for st_store qui ds `depvar' mata { // Copy data to Mata (survives any Stata data manipulation) _y_orig = st_data(., "`depvar'") _x_orig = st_data(., "`indepvar'") _res_t = st_data(., "`resid_t'") _res_f = st_data(., "`resid_f'") _TT = rows(_y_orig) // Compute original dx (first diff of x) _dx_orig = J(_TT, 1, .) for (_i = 2; _i <= _TT; _i++) _dx_orig[_i] = _x_orig[_i] - _x_orig[_i-1] // Restricted coefficients _b_t = st_matrix("`b_rest_t'") _b_f = st_matrix("`b_rest_f'") // Bootstrap storage _boot_tstat = J(`reps', 1, .) _boot_fstat = J(`reps', 1, .) // Parameters _op = `opt_p' _oq = `opt_q' _cv = `case' _si = `opt_lag' + 2 if (_si < 3) _si = 3 // Find depvar column index for st_store _depvar_idx = . _vnames = st_varindex("`depvar'") _depvar_idx = _vnames } // ========================================================================= // 4. BOOTSTRAP LOOP // Strategy: Mata generates bootstrap y, writes it to Stata via st_store, // then Stata runs regress with TS operators on the bootstrap data. // No tempfiles needed — we restore y_orig at the end. // ========================================================================= forvalues b = 1/`reps' { if mod(`b', 250) == 0 | `b' == 1 { di as txt _col(7) "Bootstrap replication `b'/`reps'..." } // ─── Mata: Recursive DGP + write to Stata ─── mata { // Joint resample with replacement _ri = ceil(uniform(_TT, 1) :* _TT) for (_i = 1; _i <= _TT; _i++) { if (_ri[_i] < 1) _ri[_i] = 1 if (_ri[_i] > _TT) _ri[_i] = _TT } _rt_b = _res_t[_ri] _rf_b = _res_f[_ri] // Initialize bootstrap series from original data _yt = _y_orig _dyt = J(_TT, 1, 0) _yf = _y_orig _dyf = J(_TT, 1, 0) for (_i = 2; _i <= _TT; _i++) { _dyt[_i] = _y_orig[_i] - _y_orig[_i-1] _dyf[_i] = _y_orig[_i] - _y_orig[_i-1] } // ──── Recursive DGP (EViews model a / a.solve) ──── for (_i = _si; _i <= _TT; _i++) { // t-test DGP: H0 pi=0, no y(-1) // Coef layout: [x(-1), dy(-1..p), dx(-1..q), det...] _val = _b_t[1] * _x_orig[_i-1] _ci = 2 for (_j = 1; _j <= _op; _j++) { if (_i-_j >= 1) _val = _val + _b_t[_ci] * _dyt[_i-_j] _ci++ } for (_k = 1; _k <= _oq; _k++) { if (_i-_k >= 2) _val = _val + _b_t[_ci] * _dx_orig[_i-_k] _ci++ } if (_cv == 3) _val = _val + _b_t[_ci] else if (_cv == 5) _val = _val + _b_t[_ci] + _b_t[_ci+1]*_i _dyt[_i] = _val + _rt_b[_i] _yt[_i] = _yt[_i-1] + _dyt[_i] // f-test DGP: H0 delta=0, no x(-1) // Coef layout: [y(-1), dy(-1..p), dx(-1..q), det...] _val = _b_f[1] * _yf[_i-1] _ci = 2 for (_j = 1; _j <= _op; _j++) { if (_i-_j >= 1) _val = _val + _b_f[_ci] * _dyf[_i-_j] _ci++ } for (_k = 1; _k <= _oq; _k++) { if (_i-_k >= 2) _val = _val + _b_f[_ci] * _dx_orig[_i-_k] _ci++ } if (_cv == 3) _val = _val + _b_f[_ci] else if (_cv == 5) _val = _val + _b_f[_ci] + _b_f[_ci+1]*_i _dyf[_i] = _val + _rf_b[_i] _yf[_i] = _yf[_i-1] + _dyf[_i] } // Write t-bootstrap y to Stata for regression st_store(., _depvar_idx, _yt) } // ─── t-test: regress on bootstrap data ─── capture { qui regress D.`depvar' `full_reg', `nocons_opt' } if _rc == 0 { local bt = _b[L.`depvar'] / _se[L.`depvar'] mata: _boot_tstat[`b'] = strtoreal(st_local("bt")) } // ─── f-test: write f-bootstrap y, then regress ─── mata: st_store(., _depvar_idx, _yf) capture { qui regress D.`depvar' `full_reg', `nocons_opt' } if _rc == 0 { capture qui test L.`indepvar' if _rc == 0 { local bf = r(F) mata: _boot_fstat[`b'] = strtoreal(st_local("bf")) } } } // ─── Restore original y ─── mata: st_store(., _depvar_idx, _y_orig) di as txt _col(7) "Bootstrap completed." di as txt "" // ========================================================================= // 5. COMPUTE CRITICAL VALUES (all in Mata) // ========================================================================= tempname tcv10 tcv05 tcv025 tcv01 fcv10 fcv05 fcv025 fcv01 nbt nbf mata { _btc = select(_boot_tstat, _boot_tstat :< .) _bfc = select(_boot_fstat, _boot_fstat :< .) _nbt = rows(_btc) _nbf = rows(_bfc) // t-test: LOWER tail (reject if observed t < cv) if (_nbt > 10) { _ts = sort(_btc, 1) st_numscalar("`tcv10'", _ts[ceil(.10*_nbt)]) st_numscalar("`tcv05'", _ts[ceil(.05*_nbt)]) st_numscalar("`tcv025'", _ts[ceil(.025*_nbt)]) st_numscalar("`tcv01'", _ts[ceil(.01*_nbt)]) } else { st_numscalar("`tcv10'", .) st_numscalar("`tcv05'", .) st_numscalar("`tcv025'", .) st_numscalar("`tcv01'", .) } // F-test: UPPER tail (reject if observed F > cv) if (_nbf > 10) { _fs = sort(_bfc, 1) st_numscalar("`fcv10'", _fs[ceil(.90*_nbf)]) st_numscalar("`fcv05'", _fs[ceil(.95*_nbf)]) st_numscalar("`fcv025'", _fs[ceil(.975*_nbf)]) st_numscalar("`fcv01'", _fs[ceil(.99*_nbf)]) } else { st_numscalar("`fcv10'", .) st_numscalar("`fcv05'", .) st_numscalar("`fcv025'", .) st_numscalar("`fcv01'", .) } st_numscalar("`nbt'", _nbt) st_numscalar("`nbf'", _nbf) } di as txt _col(5) "Valid bootstrap t-replications: " as res scalar(`nbt') di as txt _col(5) "Valid bootstrap F-replications: " as res scalar(`nbf') di as txt "" // Return all scalars return scalar t_cv10 = scalar(`tcv10') return scalar t_cv05 = scalar(`tcv05') return scalar t_cv025 = scalar(`tcv025') return scalar t_cv01 = scalar(`tcv01') return scalar f_cv10 = scalar(`fcv10') return scalar f_cv05 = scalar(`fcv05') return scalar f_cv025 = scalar(`fcv025') return scalar f_cv01 = scalar(`fcv01') return scalar B_t = scalar(`nbt') return scalar B_f = scalar(`nbf') end