*! boundeduroot_mtests v1.0.0 Merwan Roudane 07jul2026 *! GLS-based M unit-root tests for bounded processes *! Carrion-i-Silvestre & Gadea (2013), "GLS-based unit root tests for *! bounded processes" -- faithful port of the companion MATLAB code *! (estima_M_tests_bounded.m / cv_M_tests_bounded.m). *! Part of the boundeduroot library. github.com/merwanroudane program define boundeduroot_mtests, rclass version 14.0 syntax varname(ts) [if] [in] , /// Lbound(string) /// [ Ubound(string) ] /// [ Iter(integer 1000) ] /// [ MAXLag(integer -1) ] /// [ SEED(integer 16384) ] /// [ Level(cilevel) ] /// [ noGRAPH ] /// [ GNAME(string) ] * ---- bounds ------------------------------------------------------------ local lb = lower("`lbound'") local ub = lower("`ubound'") local lbval = . if !inlist("`lb'","",".","inf","-inf","none") { capture confirm number `lbound' if _rc { di as error "lbound() must be a number or . for one-sided" exit 198 } local lbval = real("`lbound'") } local ubval = . if !inlist("`ub'","",".","inf","+inf","none") { capture confirm number `ubound' if _rc { di as error "ubound() must be a number or . for one-sided" exit 198 } local ubval = real("`ubound'") } if `lbval'==. & `ubval'==. { di as error "Specify at least one finite bound in lbound() or ubound()." exit 198 } if `lbval'!=. & `ubval'!=. { if `lbval' >= `ubval' { di as error "lower bound must be strictly less than upper bound" exit 198 } } * ---- sample ------------------------------------------------------------ marksample touse markout `touse' `varlist' capture qui tsset if _rc { di as error "Data are not tsset. Use {cmd:tsset} timevar first." exit 111 } local timevar "`r(timevar)'" local panelvar "`r(panelvar)'" if "`panelvar'" != "" { di as error "Panel data are not supported; use a single time series." exit 198 } markout `touse' `timevar' qui count if `touse' local N = r(N) if `N' < 30 { di as error "Insufficient observations (need at least 30)." exit 2001 } if `maxlag' == -1 local maxlag = floor(12*(`N'/100)^0.25) if `iter' < 100 { di as error "iter() must be at least 100" exit 198 } set seed `seed' di _n as text "{bf:Bounded M unit-root tests} (Carrion-i-Silvestre & Gadea 2013)" di as text "{hline 78}" * ---- engine ------------------------------------------------------------ mata: _bmt_run("`varlist'","`touse'","`timevar'") * writes: __bmt_sta (6x3), __bmt_cv_msb __bmt_cv_mza __bmt_cv_mzt (6x4), * __bmt_c (6x2), scalars __bmt_kap_ar __bmt_kap_np __bmt_x0 tempname STA CVM CVA CVT CIN matrix `STA' = __bmt_sta matrix `CVM' = __bmt_cv_msb matrix `CVA' = __bmt_cv_mza matrix `CVT' = __bmt_cv_mzt matrix `CIN' = __bmt_c local rn "OLS_ar OLS_np GLSERS_ar GLSERS_np GLSBOUNDS_ar GLSBOUNDS_np" matrix rownames `STA' = `rn' matrix colnames `STA' = MSB MZa MZt matrix rownames `CVM' = `rn' matrix rownames `CVA' = `rn' matrix rownames `CVT' = `rn' matrix colnames `CVM' = cv1 cv2_5 cv5 cv10 matrix colnames `CVA' = cv1 cv2_5 cv5 cv10 matrix colnames `CVT' = cv1 cv2_5 cv5 cv10 matrix rownames `CIN' = `rn' matrix colnames `CIN' = c_inf c_sup * ---- configuration block ------------------------------------------------- local lbtxt = cond(`lbval'==., "-inf", string(`lbval')) local ubtxt = cond(`ubval'==., "+inf", string(`ubval')) di as text " Variable : " as result "`varlist'" di as text " Observations T : " as result `N' di as text " Bounds [b, b-bar]: " as result "[`lbtxt', `ubtxt']" di as text " X0 (first obs) : " as result %9.5f __bmt_x0 di as text " kappa-hat (GLS-BOUNDS): " as result %8.4f __bmt_kap_ar /// as text " (parametric) " as result %8.4f __bmt_kap_np as text " (nonparametric)" di as text " Simulated CVs : " as result `iter' as text " draws, folded (reflected) Brownian motion, seed(`seed')" * ---- results table ------------------------------------------------------- local plab1 "OLS demeaning" local plab2 "GLS-ERS demeaning (c-bar = -7)" local plab3 "GLS-BOUNDS demeaning (c-bar = kappa-hat)" di _n as text "{hline 78}" di as text %-26s "Detrending / LRV" _col(30) %8s "MSB" _col(44) %9s "MZa" /// _col(58) %9s "MZt" di as text "{hline 78}" forvalues p = 1/3 { di as text "{bf:`plab`p''}" forvalues s = 1/2 { local r = (`p'-1)*2 + `s' local lrvlab = cond(`s'==1, " parametric (SAR)", " nonparametric (QS)") local out "" forvalues c = 1/3 { local st = `STA'[`r',`c'] if `c'==1 matrix __cvx = `CVM' if `c'==2 matrix __cvx = `CVA' if `c'==3 matrix __cvx = `CVT' local star "" if `st' < __cvx[`r',4] local star "*" if `st' < __cvx[`r',3] local star "**" if `st' < __cvx[`r',1] local star "***" local st`c' = `st' local sr`c' "`star'" } di as text %-26s "`lrvlab'" _col(30) as result %8.4f `st1' as text %-3s "`sr1'" /// _col(44) as result %9.4f `st2' as text %-3s "`sr2'" /// _col(58) as result %9.4f `st3' as text %-3s "`sr3'" } } di as text "{hline 78}" di as text "H0: bounded unit root. All three tests reject for SMALL values of the" di as text "statistic: * p<.10 ** p<.05 *** p<.01 (bound-specific simulated CVs)." * ---- CV table ------------------------------------------------------------- local siglev = 100 - `level' local cvcol = 3 if `siglev' <= 1 local cvcol = 1 else if `siglev' <= 3 local cvcol = 2 else if `siglev' <= 5 local cvcol = 3 else local cvcol = 4 di _n as text "Simulated `siglev'% critical values (given the bounds and estimated c, c-bar):" di as text "{hline 78}" di as text %-26s "Detrending / LRV" _col(30) %8s "MSB" _col(44) %9s "MZa" /// _col(58) %9s "MZt" di as text "{hline 78}" local lab1 "OLS / SAR" local lab2 "OLS / QS" local lab3 "GLS-ERS / SAR" local lab4 "GLS-ERS / QS" local lab5 "GLS-BOUNDS / SAR" local lab6 "GLS-BOUNDS / QS" forvalues r = 1/6 { di as text %-26s " `lab`r''" _col(30) as result %8.4f `CVM'[`r',`cvcol'] /// _col(44) %9.4f `CVA'[`r',`cvcol'] _col(58) %9.4f `CVT'[`r',`cvcol'] } di as text "{hline 78}" di as text "c-hat ranges over configs: [" as result %6.3f `CIN'[1,1] as text "," /// as result %6.3f `CIN'[6,1] as text "] (lower) [" /// as result %6.3f `CIN'[1,2] as text "," as result %6.3f `CIN'[6,2] as text "] (upper)" * ---- returns --------------------------------------------------------------- return scalar N = `N' return scalar x0 = __bmt_x0 return scalar kappa_ar = __bmt_kap_ar return scalar kappa_np = __bmt_kap_np return scalar lbound = `lbval' return scalar ubound = `ubval' return scalar iter = `iter' return local depvar "`varlist'" return local timevar "`timevar'" return local cmd "boundeduroot mtests" return matrix stats = `STA', copy return matrix cv_msb = `CVM', copy return matrix cv_mza = `CVA', copy return matrix cv_mzt = `CVT', copy return matrix cpars = `CIN', copy * ---- graph ------------------------------------------------------------------- if "`graph'" != "nograph" { if "`gname'" == "" local gname bmtests _bmt_plot , timevar(`timevar') depvar(`varlist') touse(`touse') /// lbound(`lbval') ubound(`ubval') gname(`gname') } capture matrix drop __bmt_sta __bmt_cv_msb __bmt_cv_mza __bmt_cv_mzt __bmt_c __cvx capture scalar drop __bmt_kap_ar __bmt_kap_np __bmt_x0 end *============================================================================== * Journal figure: series + bounds, and statistic-vs-CV dot chart *============================================================================== program define _bmt_plot version 14.0 syntax , timevar(string) depvar(string) touse(string) /// [ lbound(string) ubound(string) gname(string) ] local yl "" if !inlist("`lbound'",".","") /// local yl "`yl' yline(`lbound', lpattern(dash) lcolor(red) lwidth(medthin))" if !inlist("`ubound'",".","") /// local yl "`yl' yline(`ubound', lpattern(dash) lcolor(red) lwidth(medthin))" twoway (line `depvar' `timevar' if `touse', lcolor(navy) lwidth(medthin)), /// `yl' /// title("Bounded series", size(medium)) /// subtitle("dashed red = bounds", size(small)) /// ytitle("`depvar'") xtitle("`timevar'") /// graphregion(color(white)) plotregion(color(white)) /// name(`gname'_series, replace) nodraw * dot chart: statistic vs 5% critical value across the 6 configurations preserve qui { clear set obs 6 gen cfg = _n gen msb = . gen mza = . gen mzt = . gen cvmsb = . gen cvmza = . gen cvmzt = . forvalues r = 1/6 { replace msb = __bmt_sta[`r',1] in `r' replace mza = __bmt_sta[`r',2] in `r' replace mzt = __bmt_sta[`r',3] in `r' replace cvmsb = __bmt_cv_msb[`r',3] in `r' replace cvmza = __bmt_cv_mza[`r',3] in `r' replace cvmzt = __bmt_cv_mzt[`r',3] in `r' } label define __bmtc 1 "OLS/SAR" 2 "OLS/QS" 3 "ERS/SAR" 4 "ERS/QS" /// 5 "BND/SAR" 6 "BND/QS", replace label values cfg __bmtc } twoway (scatter mzt cfg, mcolor(navy) msymbol(circle)) /// (scatter cvmzt cfg, mcolor(red) msymbol(diamond)), /// legend(order(1 "MZt statistic" 2 "5% simulated CV") rows(1) size(small)) /// title("MZt vs bound-specific critical value", size(medium)) /// xlabel(1(1)6, valuelabel angle(45) labsize(small)) /// xtitle("") ytitle("MZt") /// graphregion(color(white)) plotregion(color(white)) /// name(`gname'_cv, replace) nodraw restore capture graph combine `gname'_series `gname'_cv, /// title("boundeduroot mtests: Carrion-i-Silvestre & Gadea (2013)", size(medsmall)) /// note("Below its CV the statistic rejects the bounded unit-root null.", size(vsmall)) /// graphregion(color(white)) name(`gname', replace) if _rc { capture graph display `gname'_series } end *============================================================================== * M A T A E N G I N E (_bmt_*) *============================================================================== version 14.0 mata: // --------------------------------------------------------------------------- // Skorohod folding of a path into [binf, bsup] (rbm.m) // A missing bound (.) disables that side. // --------------------------------------------------------------------------- real colvector _bmt_rbm(real colvector x0, real scalar binf, real scalar bsup) { real colvector x real scalar it, bad x = x0 it = 0 bad = 1 while (bad & it < 10000) { bad = 0 if (bsup < .) { if (max(x) > bsup) { x = bsup :- abs(x :- bsup) bad = 1 } } if (binf < .) { if (min(x) < binf) { x = binf :+ abs(x :- binf) bad = 1 } } it = it + 1 } return(x) } // --------------------------------------------------------------------------- // GLS (quasi-difference) demeaning with noncentrality cbar (det_gls) // --------------------------------------------------------------------------- real colvector _bmt_detgls(real colvector y, real scalar cbar) { real scalar nt, abar, bhat real colvector ya, za, z nt = rows(y) z = J(nt,1,1) abar = 1 + cbar/nt ya = y za = z ya[|2 \ nt|] = y[|2 \ nt|] :- abar*y[|1 \ nt-1|] za[|2 \ nt|] = z[|2 \ nt|] :- abar*z[|1 \ nt-1|] bhat = quadcross(za,ya)/quadcross(za,za) return(y :- z*bhat) } // --------------------------------------------------------------------------- // Ng-Perron MAIC lag selection on a detrended series (s2ar) // --------------------------------------------------------------------------- real scalar _bmt_maic(real colvector d, real scalar kmax) { real scalar T, nef, k, i, r, t, sumy, s2, tau, mic, best, bestv real colvector dep, bb, e real matrix Rg, Xk T = rows(d) nef = T - kmax - 1 if (nef < 5) return(0) dep = J(nef,1,0) Rg = J(nef, kmax+1, 0) for (r=1; r<=nef; r++) { t = kmax + 1 + r dep[r] = d[t] - d[t-1] Rg[r,1] = d[t-1] for (i=1; i<=kmax; i++) { Rg[r,1+i] = d[t-i] - d[t-i-1] } } sumy = quadcross(Rg[.,1], Rg[.,1]) best = 0 bestv = . for (k=0; k<=kmax; k++) { Xk = Rg[|1,1 \ nef,k+1|] bb = invsym(quadcross(Xk,Xk))*quadcross(Xk,dep) e = dep - Xk*bb s2 = quadcross(e,e)/nef tau = (bb[1]^2 * sumy)/s2 mic = ln(s2) + 2*(k+tau)/nef if (mic < bestv | k==0) { bestv = mic best = k } } return(best) } // --------------------------------------------------------------------------- // AR spectral long-run variance from ADF regression with k lags (adfp) // on an already-detrended series; s2vec = s2e/(1-sum b)^2. // --------------------------------------------------------------------------- real scalar _bmt_sar(real colvector d, real scalar k) { real scalar T, nobs, r, t, i, s2, sumb real colvector dep, b, e real matrix Rg T = rows(d) nobs = T - k - 1 dep = J(nobs,1,0) Rg = J(nobs, k+1, 0) for (r=1; r<=nobs; r++) { t = k + 1 + r dep[r] = d[t] - d[t-1] Rg[r,1] = d[t-1] for (i=1; i<=k; i++) { Rg[r,1+i] = d[t-i] - d[t-i-1] } } b = invsym(quadcross(Rg,Rg))*quadcross(Rg,dep) e = dep - Rg*b s2 = quadcross(e,e)/nobs sumb = 0 if (k >= 1) sumb = sum(b[|2 \ k+1|]) return(s2/((1-sumb)^2)) } // --------------------------------------------------------------------------- // Parametric long-run variance (estima_lrv_ar): // detm==0 : OLS demeaning (Perron-Qu 2007) // detm==1 : GLS demeaning at cbar (Ng-Perron 2001) // --------------------------------------------------------------------------- real scalar _bmt_lrv_ar(real colvector y, real scalar cbar, real scalar detm, real scalar kmax) { real colvector d real scalar kopt if (detm == 0) { d = y :- mean(y) } else { d = _bmt_detgls(y, cbar) } kopt = _bmt_maic(d, kmax) return(_bmt_sar(d, kopt)) } // --------------------------------------------------------------------------- // Nonparametric long-run variance (estima_lrv_np), QS kernel with the // Newey-West (1994) automatic bandwidth. // --------------------------------------------------------------------------- real scalar _bmt_lrv_np(real colvector y, real scalar cbar, real scalar detm) { real colvector d, dy, res, acov real scalar T, nr, beta, n, j, s0, s2w, gam, mbw, xx, kw, lrv T = rows(y) if (detm == 0) { d = y :- mean(y) } else { d = _bmt_detgls(y, cbar) } dy = d[|2 \ T|] - d[|1 \ T-1|] beta = quadcross(d[|1 \ T-1|], dy)/quadcross(d[|1 \ T-1|], d[|1 \ T-1|]) res = dy - d[|1 \ T-1|]*beta nr = rows(res) acov = J(nr,1,0) for (j=0; j<=nr-1; j++) { acov[j+1] = quadcross(res[|j+1 \ nr|], res[|1 \ nr-j|])/nr } n = floor(4*(T/100)^(2/25)) if (n > nr-1) n = nr-1 mbw = 0 if (n > 0) { s0 = acov[1] s2w = 0 for (j=1; j<=n; j++) { s0 = s0 + 2*acov[j+1] s2w = s2w + 2*(j^2)*acov[j+1] } if (s0 != 0) { gam = 1.3221*((s2w/s0)^2)^(1/5) mbw = gam*T^(1/5) if (mbw > T) mbw = T } } lrv = acov[1] if (mbw > 0) { for (j=1; j<=nr-1; j++) { xx = j/mbw kw = (25/(12*pi()^2*xx^2))*(sin(1.2*pi()*xx)/(1.2*pi()*xx) - cos(1.2*pi()*xx)) lrv = lrv + 2*acov[j+1]*kw } } if (lrv <= 0) lrv = acov[1] return(lrv) } // --------------------------------------------------------------------------- // Bilinear interpolation of kappa(c_low, c_up) on the CSG (2013) grid // (read_values_kappa, 417 points). Values clamped to the grid hull. // --------------------------------------------------------------------------- real matrix _bmt_kappa_table() { string scalar kstr real colvector v real matrix K real scalar i kstr = "" kstr = kstr + "-1.0500 0 -7.2188 -1.0500 0.1000 -8.5000 -1.0500 0.1500 -11.0625 -1.0500 0.2000 -15.3750 -1.0500 0.2500 -16.5625 -1.0500 0.3000 -14.4160 -1.0500 0.3500 -12.3395 -1.0500 0.4000 -10.5712 " kstr = kstr + "-1.0500 0.4500 -9.4609 -1.0500 0.5000 -8.6716 -1.0500 0.5500 -8.0278 -1.0500 0.6000 -7.6946 -1.0500 0.6500 -7.4791 -1.0500 0.7000 -7.3560 -1.0500 0.7500 -7.2637 -1.0500 0.8000 -7.2175 " kstr = kstr + "-1.0500 0.8500 -7.2175 -1.0500 0.9000 -7.2175 -1.0500 0.9500 -7.2175 -1.0500 1.0000 -7.2175 -1.0500 1.0500 -7.2175 -1.0000 0 -7.2175 -1.0000 0.1000 -8.4983 -1.0000 0.1500 -11.0608 " kstr = kstr + "-1.0000 0.2000 -15.3733 -1.0000 0.2500 -16.5921 -1.0000 0.3000 -14.4154 -1.0000 0.3500 -12.3493 -1.0000 0.4000 -10.5581 -1.0000 0.4500 -9.4623 -1.0000 0.5000 -8.6589 -1.0000 0.5500 -8.0283 " kstr = kstr + "-1.0000 0.6000 -7.6859 -1.0000 0.6500 -7.4804 -1.0000 0.7000 -7.3588 -1.0000 0.7500 -7.2511 -1.0000 0.8000 -7.2242 -1.0000 0.8500 -7.2126 -1.0000 0.9000 -7.2126 -1.0000 0.9500 -7.2126 " kstr = kstr + "-1.0000 1.0000 -7.2126 -1.0000 1.0500 -7.2126 -0.9500 0 -7.2126 -0.9500 0.1000 -8.5016 -0.9500 0.1500 -11.0328 -0.9500 0.2000 -15.4078 -0.9500 0.2500 -16.5953 -0.9500 0.3000 -14.4215 " kstr = kstr + "-0.9500 0.3500 -12.3303 -0.9500 0.4000 -10.5633 -0.9500 0.4500 -9.4622 -0.9500 0.5000 -8.6635 -0.9500 0.5500 -8.0327 -0.9500 0.6000 -7.6904 -0.9500 0.6500 -7.4849 -0.9500 0.7000 -7.3479 " kstr = kstr + "-0.9500 0.7500 -7.2554 -0.9500 0.8000 -7.2094 -0.9500 0.8500 -7.2094 -0.9500 0.9000 -7.2094 -0.9500 0.9500 -7.2094 -0.9500 1.0000 -7.2094 -0.9500 1.0500 -7.2094 -0.9000 0 -7.2094 " kstr = kstr + "-0.9000 0.1000 -8.5058 -0.9000 0.1500 -11.0371 -0.9000 0.2000 -15.4121 -0.9000 0.2500 -16.5996 -0.9000 0.3000 -14.4258 -0.9000 0.3500 -12.3346 -0.9000 0.4000 -10.5675 -0.9000 0.4500 -9.4573 " kstr = kstr + "-0.9000 0.5000 -8.6680 -0.9000 0.5500 -8.0317 -0.9000 0.6000 -7.6908 -0.9000 0.6500 -7.4848 -0.9000 0.7000 -7.3478 -0.9000 0.7500 -7.2553 -0.9000 0.8000 -7.2094 -0.9000 0.8500 -7.2094 " kstr = kstr + "-0.9000 0.9000 -7.2094 -0.9000 0.9500 -7.2094 -0.9000 1.0000 -7.2094 -0.9000 1.0500 -7.2094 -0.8500 0 -7.2094 -0.8500 0.1000 -8.5058 -0.8500 0.1500 -11.0370 -0.8500 0.2000 -15.4120 " kstr = kstr + "-0.8500 0.2500 -16.5995 -0.8500 0.3000 -14.4257 -0.8500 0.3500 -12.3345 -0.8500 0.4000 -10.5675 -0.8500 0.4500 -9.4572 -0.8500 0.5000 -8.6679 -0.8500 0.5500 -8.0316 -0.8500 0.6000 -7.6907 " kstr = kstr + "-0.8500 0.6500 -7.4848 -0.8500 0.7000 -7.3477 -0.8500 0.7500 -7.2553 -0.8500 0.8000 -7.2093 -0.8500 0.8500 -7.2093 -0.8500 0.9000 -7.2093 -0.8500 0.9500 -7.2093 -0.8500 1.0000 -7.2093 " kstr = kstr + "-0.8500 1.0500 -7.2093 -0.8000 0 -7.2093 -0.8000 0.1000 -8.5057 -0.8000 0.1500 -11.0370 -0.8000 0.2000 -15.4120 -0.8000 0.2500 -16.5995 -0.8000 0.3000 -14.4256 -0.8000 0.3500 -12.3344 " kstr = kstr + "-0.8000 0.4000 -10.5674 -0.8000 0.4500 -9.4572 -0.8000 0.5000 -8.6679 -0.8000 0.5500 -8.0316 -0.8000 0.6000 -7.6907 -0.8000 0.6500 -7.4847 -0.8000 0.7000 -7.3477 -0.8000 0.7500 -7.2552 " kstr = kstr + "-0.8000 0.8000 -7.2092 -0.8000 0.8500 -7.2092 -0.8000 0.9000 -7.2092 -0.8000 0.9500 -7.2092 -0.8000 1.0000 -7.2092 -0.8000 1.0500 -7.2092 -0.7500 0 -7.2556 -0.7500 0.1000 -8.5056 " kstr = kstr + "-0.7500 0.1500 -11.0369 -0.7500 0.2000 -15.4119 -0.7500 0.2500 -16.5994 -0.7500 0.3000 -14.4256 -0.7500 0.3500 -12.3344 -0.7500 0.4000 -10.5674 -0.7500 0.4500 -9.4571 -0.7500 0.5000 -8.6678 " kstr = kstr + "-0.7500 0.5500 -8.0315 -0.7500 0.6000 -7.6906 -0.7500 0.6500 -7.4846 -0.7500 0.7000 -7.3630 -0.7500 0.7500 -7.2553 -0.7500 0.8000 -7.2553 -0.7500 0.8500 -7.2094 -0.7500 0.9000 -7.2094 " kstr = kstr + "-0.7500 0.9500 -7.2094 -0.7500 1.0000 -7.2094 -0.7500 1.0500 -7.2094 -0.7000 0 -7.3808 -0.7000 0.1000 -8.5058 -0.7000 0.1500 -11.0370 -0.7000 0.2000 -15.4120 -0.7000 0.2500 -16.5995 " kstr = kstr + "-0.7000 0.3000 -14.4257 -0.7000 0.3500 -12.3345 -0.7000 0.4000 -10.5675 -0.7000 0.4500 -9.4572 -0.7000 0.5000 -8.6680 -0.7000 0.5500 -8.0620 -0.7000 0.6000 -7.7286 -0.7000 0.6500 -7.5212 " kstr = kstr + "-0.7000 0.7000 -7.3706 -0.7000 0.7500 -7.3099 -0.7000 0.8000 -7.2346 -0.7000 0.8500 -7.2346 -0.7000 0.9000 -7.2346 -0.7000 0.9500 -7.2346 -0.7000 1.0000 -7.2346 -0.7000 1.0500 -7.2346 " kstr = kstr + "-0.6500 0 -7.6711 -0.6500 0.1000 -8.5305 -0.6500 0.1500 -11.0617 -0.6500 0.2000 -15.3742 -0.6500 0.2500 -16.5617 -0.6500 0.3000 -14.4152 -0.6500 0.3500 -12.3387 -0.6500 0.4000 -10.5704 " kstr = kstr + "-0.6500 0.4500 -9.4785 -0.6500 0.5000 -8.7239 -0.6500 0.5500 -8.1484 -0.6500 0.6000 -7.8333 -0.6500 0.6500 -7.6670 -0.6500 0.7000 -7.5288 -0.6500 0.7500 -7.4826 -0.6500 0.8000 -7.4518 " kstr = kstr + "-0.6500 0.8500 -7.4441 -0.6500 0.9000 -7.4441 -0.6500 0.9500 -7.4441 -0.6500 1.0000 -7.4441 -0.6500 1.0500 -7.4441 -0.6000 0 -8.0456 -0.6000 0.1000 -8.7331 -0.6000 0.1500 -11.0456 " kstr = kstr + "-0.6000 0.2000 -15.3893 -0.6000 0.2500 -16.5768 -0.6000 0.3000 -14.4079 -0.6000 0.3500 -12.3384 -0.6000 0.4000 -10.7014 -0.6000 0.4500 -9.5872 -0.6000 0.5000 -8.8638 -0.6000 0.5500 -8.3485 " kstr = kstr + "-0.6000 0.6000 -8.0456 -0.6000 0.6500 -7.8822 -0.6000 0.7000 -7.8060 -0.6000 0.7500 -7.7752 -0.6000 0.8000 -7.7752 -0.6000 0.8500 -7.7598 -0.6000 0.9000 -7.7598 -0.6000 0.9500 -7.7598 " kstr = kstr + "-0.6000 1.0000 -7.7598 -0.6000 1.0500 -7.7598 -0.5500 0 -8.6502 -0.5500 0.1000 -9.0877 -0.5500 0.1500 -11.2439 -0.5500 0.2000 -15.4627 -0.5500 0.2500 -16.5877 -0.5500 0.3000 -14.4598 " kstr = kstr + "-0.5500 0.3500 -12.4291 -0.5500 0.4000 -10.7445 -0.5500 0.4500 -9.8055 -0.5500 0.5000 -9.1283 -0.5500 0.5500 -8.6439 -0.5500 0.6000 -8.4318 -0.5500 0.6500 -8.3125 -0.5500 0.7000 -8.2391 " kstr = kstr + "-0.5500 0.7500 -8.1936 -0.5500 0.8000 -8.1782 -0.5500 0.8500 -8.1782 -0.5500 0.9000 -8.1782 -0.5500 0.9500 -8.1782 -0.5500 1.0000 -8.1782 -0.5500 1.0500 -8.1782 -0.5000 0 -9.4436 " kstr = kstr + "-0.5000 0.1000 -9.5373 -0.5000 0.1500 -11.8811 -0.5000 0.2000 -15.6936 -0.5000 0.2500 -16.6311 -0.5000 0.3000 -14.5031 -0.5000 0.3500 -12.5417 -0.5000 0.4000 -11.0461 -0.5000 0.4500 -10.1076 " kstr = kstr + "-0.5000 0.5000 -9.4302 -0.5000 0.5500 -9.0369 -0.5000 0.6000 -8.8607 -0.5000 0.6500 -8.7467 -0.5000 0.7000 -8.7467 -0.5000 0.7500 -8.7171 -0.5000 0.8000 -8.7171 -0.5000 0.8500 -8.7171 " kstr = kstr + "-0.5000 0.9000 -8.7171 -0.5000 0.9500 -8.7171 -0.5000 1.0000 -8.7171 -0.5000 1.0500 -8.7171 -0.4500 0 -10.5904 -0.4500 0.1000 -10.4078 -0.4500 0.1500 -12.8131 -0.4500 0.2000 -16.8756 " kstr = kstr + "-0.4500 0.2500 -16.9693 -0.4500 0.3000 -14.7579 -0.4500 0.3500 -12.9347 -0.4500 0.4000 -11.4871 -0.4500 0.4500 -10.6041 -0.4500 0.5000 -10.1008 -0.4500 0.5500 -9.7504 -0.4500 0.6000 -9.6619 " kstr = kstr + "-0.4500 0.6500 -9.6174 -0.4500 0.7000 -9.5714 -0.4500 0.7500 -9.5560 -0.4500 0.8000 -9.5099 -0.4500 0.8500 -9.5099 -0.4500 0.9000 -9.5099 -0.4500 0.9500 -9.5099 -0.4500 1.0000 -9.5099 " kstr = kstr + "-0.4500 1.0500 -9.5099 -0.4000 0 -12.2594 -0.4000 0.1000 -11.6227 -0.4000 0.1500 -14.4313 -0.4000 0.2000 -18.4000 -0.4000 0.2500 -17.6725 -0.4000 0.3000 -15.3953 -0.4000 0.3500 -13.5595 " kstr = kstr + "-0.4000 0.4000 -12.2651 -0.4000 0.4500 -11.4767 -0.4000 0.5000 -11.0677 -0.4000 0.5500 -10.8063 -0.4000 0.6000 -10.7601 -0.4000 0.6500 -10.7140 -0.4000 0.7000 -10.6986 -0.4000 0.7500 -10.6986 " kstr = kstr + "-0.4000 0.8000 -10.6986 -0.4000 0.8500 -10.6986 -0.4000 0.9000 -10.6986 -0.4000 0.9500 -10.6986 -0.4000 1.0000 -10.6986 -0.4000 1.0500 -10.6986 -0.3500 0 -14.5265 -0.3500 0.1000 -13.7023 " kstr = kstr + "-0.3500 0.1500 -17.0421 -0.3500 0.2000 -20.6359 -0.3500 0.2500 -18.8292 -0.3500 0.3000 -16.5165 -0.3500 0.3500 -14.5341 -0.3500 0.4000 -13.3634 -0.3500 0.4500 -12.7256 -0.3500 0.5000 -12.3300 " kstr = kstr + "-0.3500 0.5500 -12.1865 -0.3500 0.6000 -12.0963 -0.3500 0.6500 -12.0963 -0.3500 0.7000 -12.0660 -0.3500 0.7500 -12.0660 -0.3500 0.8000 -12.0660 -0.3500 0.8500 -12.0660 -0.3500 0.9000 -12.0660 " kstr = kstr + "-0.3500 0.9500 -12.0660 -0.3500 1.0000 -12.0660 -0.3500 1.0500 -12.0660 -0.3000 0 -17.9401 -0.3000 0.1000 -16.7467 -0.3000 0.1500 -21.6862 -0.3000 0.2000 -23.3112 -0.3000 0.2500 -20.8072 " kstr = kstr + "-0.3000 0.3000 -17.9358 -0.3000 0.3500 -16.2863 -0.3000 0.4000 -15.2209 -0.3000 0.4500 -14.6005 -0.3000 0.5000 -14.2680 -0.3000 0.5500 -14.2376 -0.3000 0.6000 -14.2376 -0.3000 0.6500 -14.2376 " kstr = kstr + "-0.3000 0.7000 -14.2376 -0.3000 0.7500 -14.2376 -0.3000 0.8000 -14.2376 -0.3000 0.8500 -14.2376 -0.3000 0.9000 -14.2376 -0.3000 0.9500 -14.2376 -0.3000 1.0000 -14.2376 -0.3000 1.0500 -14.2376 " kstr = kstr + "-0.2500 0 -23.2055 -0.2500 0.1000 -22.1723 -0.2500 0.1500 -28.4555 -0.2500 0.2000 -27.5023 -0.2500 0.2500 -23.1836 -0.2500 0.3000 -20.5154 -0.2500 0.3500 -18.8388 -0.2500 0.4000 -17.7081 " kstr = kstr + "-0.2500 0.4500 -17.0341 -0.2500 0.5000 -16.8697 -0.2500 0.5500 -16.7255 -0.2500 0.6000 -16.7107 -0.2500 0.6500 -16.7107 -0.2500 0.7000 -16.7107 -0.2500 0.7500 -16.7107 -0.2500 0.8000 -16.7107 " kstr = kstr + "-0.2500 0.8500 -16.7107 -0.2500 0.9000 -16.7107 -0.2500 0.9500 -16.7107 -0.2500 1.0000 -16.7107 -0.2500 1.0500 -16.7107 -0.2000 0 -32.4755 -0.2000 0.1000 -32.6005 -0.2000 0.1500 -38.4755 " kstr = kstr + "-0.2000 0.2000 -32.4616 -0.2000 0.2500 -27.3497 -0.2000 0.3000 -23.4455 -0.2000 0.3500 -20.4055 -0.2000 0.4000 -18.2375 -0.2000 0.4500 -16.6459 -0.2000 0.5000 -15.9442 -0.2000 0.5500 -15.5319 " kstr = kstr + "-0.2000 0.6000 -15.4731 -0.2000 0.6500 -15.4731 -0.2000 0.7000 -15.4731 -0.2000 0.7500 -15.4731 -0.2000 0.8000 -15.4731 -0.2000 0.8500 -15.4731 -0.2000 0.9000 -15.4731 -0.2000 0.9500 -15.4731 " kstr = kstr + "-0.2000 1.0000 -15.4731 -0.2000 1.0500 -15.4731 -0.1500 0 -51.2194 -0.1500 0.1000 -55.9069 -0.1500 0.1500 -51.2265 -0.1500 0.2000 -38.2680 -0.1500 0.2500 -28.3652 -0.1500 0.3000 -21.4902 " kstr = kstr + "-0.1500 0.3500 -17.2152 -0.1500 0.4000 -14.7971 -0.1500 0.4500 -12.8630 -0.1500 0.5000 -12.0242 -0.1500 0.5500 -11.4961 -0.1500 0.6000 -11.2067 -0.1500 0.6500 -11.0883 -0.1500 0.7000 -11.0883 " kstr = kstr + "-0.1500 0.7500 -11.0883 -0.1500 0.8000 -11.0883 -0.1500 0.8500 -11.0883 -0.1500 0.9000 -11.0883 -0.1500 0.9500 -11.0883 -0.1500 1.0000 -11.0883 -0.1500 1.0500 -11.0883 -0.1000 0 -99.5817 " kstr = kstr + "-0.1000 0.1000 -99.5817 -0.1000 0.1500 -55.4845 -0.1000 0.2000 -31.9151 -0.1000 0.2500 -22.0123 -0.1000 0.3000 -16.4706 -0.1000 0.3500 -13.5372 -0.1000 0.4000 -11.6924 -0.1000 0.4500 -10.5616 " kstr = kstr + "-0.1000 0.5000 -9.4556 -0.1000 0.5500 -8.8475 -0.1000 0.6000 -8.6355 -0.1000 0.6500 -8.4868 -0.1000 0.7000 -8.4423 -0.1000 0.7500 -8.3655 -0.1000 0.8000 -8.3348 -0.1000 0.8500 -8.3348 " kstr = kstr + "-0.1000 0.9000 -8.3348 " v = strtoreal(tokens(kstr))' K = J(rows(v)/3, 3, .) for (i=1; i<=rows(K); i++) { K[i,1] = v[3*i-2] K[i,2] = v[3*i-1] K[i,3] = v[3*i] } return(K) } real scalar _bmt_kappa(real scalar cinf0, real scalar csup0) { real matrix K real colvector clg, cug real scalar cinf, csup, il, iu, i, w1, w2 real scalar k11, k12, k21, k22, kl, kh K = _bmt_kappa_table() // grids: c_low ascending -1.05..-0.10 step .05; c_up {0,.10,.15,...,1.05} clg = range(-1.05, -0.10, 0.05) cug = 0 \ range(0.10, 1.05, 0.05) cinf = cinf0 csup = csup0 if (cinf >= .) cinf = -1.05 if (csup >= .) csup = 1.05 if (cinf < -1.05) cinf = -1.05 if (cinf > -0.10) cinf = -0.10 if (csup < 0.00) csup = 0.00 if (csup > 1.05) csup = 1.05 il = 1 for (i=1; i<=rows(clg)-1; i++) { if (cinf >= clg[i]) il = i } iu = 1 for (i=1; i<=rows(cug)-1; i++) { if (csup >= cug[i]) iu = i } k11 = _bmt_klook(K, clg[il], cug[iu]) k12 = _bmt_klook(K, clg[il], cug[iu+1]) k21 = _bmt_klook(K, clg[il+1], cug[iu]) k22 = _bmt_klook(K, clg[il+1], cug[iu+1]) w1 = 0 if (clg[il+1] > clg[il]) w1 = (cinf - clg[il])/(clg[il+1] - clg[il]) w2 = 0 if (cug[iu+1] > cug[iu]) w2 = (csup - cug[iu])/(cug[iu+1] - cug[iu]) kl = k11 + w2*(k12 - k11) kh = k21 + w2*(k22 - k21) return(kl + w1*(kh - kl)) } real scalar _bmt_klook(real matrix K, real scalar cl, real scalar cu) { real scalar i, best, bestd, d best = . bestd = . for (i=1; i<=rows(K); i++) { if (abs(K[i,1]-cl) < 1e-6) { d = abs(K[i,2]-cu) if (d < bestd | bestd >= .) { bestd = d best = K[i,3] } } } return(best) } // --------------------------------------------------------------------------- // M statistics (2013 convention: no X0^2 term in the MZa numerator) // returns (MSB, MZa, MZt) // --------------------------------------------------------------------------- real rowvector _bmt_mstats(real colvector d, real scalar lrv) { real scalar T, ss, msb, mza T = rows(d) ss = quadcross(d[|1 \ T-1|], d[|1 \ T-1|]) msb = sqrt(ss/(lrv*T^2)) mza = (d[T]^2/T - lrv)/(2*ss/T^2) return((msb, mza, mza*msb)) } real scalar _bmt_quantile(real colvector x, real scalar p) { real colvector s real scalar n, h, fl s = sort(x,1) n = rows(s) h = (n-1)*p + 1 fl = floor(h) if (fl >= n) return(s[n]) if (fl < 1) return(s[1]) return(s[fl] + (h-fl)*(s[fl+1]-s[fl])) } // --------------------------------------------------------------------------- // Driver // --------------------------------------------------------------------------- void _bmt_run(string scalar yvar, string scalar touse, string scalar tvar) { real matrix X, STA, CVM, CVA, CVT, C real colvector y, d, x, ysim, dsim real rowvector mm, kapcfg, detcfg, npcfg real scalar T, kmax, iter, lb, ub, x0, i, j, c real scalar lrv, lrv0, cinf0, csup0, kap0, kap_ar, kap_np, kk real scalar cinf, csup, lsim real matrix S string scalar rs lb = strtoreal(st_local("lbval")) ub = strtoreal(st_local("ubval")) kmax = strtoreal(st_local("maxlag")) iter = strtoreal(st_local("iter")) X = st_data(., (yvar, tvar), touse) X = sort(X, 2) y = X[.,1] T = rows(y) x0 = y[1] STA = J(6,3,.) C = J(6,2,.) CVM = J(6,4,.) CVA = J(6,4,.) CVT = J(6,4,.) // ------- observed statistics, config by config (estima_M_tests_bounded) // 1: OLS + SAR d = y :- mean(y) lrv = _bmt_lrv_ar(y, ., 0, kmax) STA[1,.] = _bmt_mstats(d, lrv) C[1,1] = _bmt_cpar(lb, x0, lrv, T) C[1,2] = _bmt_cpar(ub, x0, lrv, T) // 2: OLS + NP lrv = _bmt_lrv_np(y, ., 0) STA[2,.] = _bmt_mstats(d, lrv) C[2,1] = _bmt_cpar(lb, x0, lrv, T) C[2,2] = _bmt_cpar(ub, x0, lrv, T) // 3: GLS-ERS + SAR d = _bmt_detgls(y, -7) lrv = _bmt_lrv_ar(y, -7, 1, kmax) STA[3,.] = _bmt_mstats(d, lrv) C[3,1] = _bmt_cpar(lb, x0, lrv, T) C[3,2] = _bmt_cpar(ub, x0, lrv, T) // 4: GLS-ERS + NP lrv = _bmt_lrv_np(y, -7, 1) STA[4,.] = _bmt_mstats(d, lrv) C[4,1] = _bmt_cpar(lb, x0, lrv, T) C[4,2] = _bmt_cpar(ub, x0, lrv, T) // 5: GLS-BOUNDS + SAR (two-stage kappa, stats use stage-1 lrv) lrv0 = _bmt_lrv_ar(y, ., 0, kmax) cinf0 = _bmt_cpar(lb, x0, lrv0, T) csup0 = _bmt_cpar(ub, x0, lrv0, T) kap0 = _bmt_kappa(cinf0, csup0) lrv = _bmt_lrv_ar(y, kap0, 1, kmax) C[5,1] = _bmt_cpar(lb, x0, lrv, T) C[5,2] = _bmt_cpar(ub, x0, lrv, T) kap_ar = _bmt_kappa(C[5,1], C[5,2]) d = _bmt_detgls(y, kap_ar) STA[5,.] = _bmt_mstats(d, lrv) // 6: GLS-BOUNDS + NP lrv0 = _bmt_lrv_np(y, ., 0) cinf0 = _bmt_cpar(lb, x0, lrv0, T) csup0 = _bmt_cpar(ub, x0, lrv0, T) kap0 = _bmt_kappa(cinf0, csup0) lrv = _bmt_lrv_np(y, kap0, 1) C[6,1] = _bmt_cpar(lb, x0, lrv, T) C[6,2] = _bmt_cpar(ub, x0, lrv, T) kap_np = _bmt_kappa(C[6,1], C[6,2]) d = _bmt_detgls(y, kap_np) STA[6,.] = _bmt_mstats(d, lrv) // ------- simulated critical values (cv_M_tests_bounded) ---------------- // config -> (detrending kappa or OLS, parametric/np) detcfg = (0, 0, 1, 1, 1, 1) // 0 OLS, 1 GLS kapcfg = (., ., -7, -7, _bmt_kappa(C[5,1],C[5,2]), _bmt_kappa(C[6,1],C[6,2])) npcfg = (0, 1, 0, 1, 0, 1) // 0 SAR, 1 NP rs = rseed() for (c=1; c<=6; c++) { rseed(rs) // common random-walk bank across configs cinf = C[c,1] csup = C[c,2] S = J(iter,3,.) for (j=1; j<=iter; j++) { x = J(T,1,0) x[|2 \ T|] = rnormal(T-1,1,0,1) x = quadrunningsum(x) ysim = _bmt_rbm(x, _bmt_scb(cinf,T), _bmt_scb(csup,T)) if (detcfg[c] == 0) { dsim = ysim :- mean(ysim) } else { dsim = _bmt_detgls(ysim, kapcfg[c]) } if (npcfg[c] == 0) { lsim = _bmt_lrv_ar(ysim, kapcfg[c], detcfg[c], kmax) } else { lsim = _bmt_lrv_np(ysim, kapcfg[c], detcfg[c]) } S[j,.] = _bmt_mstats(dsim, lsim) } CVM[c,.] = (_bmt_quantile(S[.,1],.01), _bmt_quantile(S[.,1],.025), _bmt_quantile(S[.,1],.05), _bmt_quantile(S[.,1],.10)) CVA[c,.] = (_bmt_quantile(S[.,2],.01), _bmt_quantile(S[.,2],.025), _bmt_quantile(S[.,2],.05), _bmt_quantile(S[.,2],.10)) CVT[c,.] = (_bmt_quantile(S[.,3],.01), _bmt_quantile(S[.,3],.025), _bmt_quantile(S[.,3],.05), _bmt_quantile(S[.,3],.10)) } st_matrix("__bmt_sta", STA) st_matrix("__bmt_cv_msb", CVM) st_matrix("__bmt_cv_mza", CVA) st_matrix("__bmt_cv_mzt", CVT) st_matrix("__bmt_c", C) st_numscalar("__bmt_kap_ar", kap_ar) st_numscalar("__bmt_kap_np", kap_np) st_numscalar("__bmt_x0", x0) } // bound parameter c = (b - X0)/sqrt(lrv*T); missing bound stays missing real scalar _bmt_cpar(real scalar b, real scalar x0, real scalar lrv, real scalar T) { if (b >= .) return(.) return((b - x0)/sqrt(lrv*T)) } // scaled bound c*sqrt(T); missing stays missing (disables folding side) real scalar _bmt_scb(real scalar cv, real scalar T) { if (cv >= .) return(.) return(cv*sqrt(T)) } end