*! _qnardl_bounds v0.3.0 27may2026 *! Author: Dr Merwan Roudane *! PSS bounds testing at each tau, per Bertsatos, Sakellaris & Tsionas (2022). *! *! Three tests per tau (with case-dependent critical values): *! F_yx : joint H0 phi_y(tau) = 0 AND phi+_j(tau) = phi-_j(tau) = delta_j = 0 *! F_x : joint H0 phi+_j(tau) = phi-_j(tau) = delta_j = 0 (excluding phi_y) *! t_y : H0 phi_y(tau) = 0 (one-sided, expected negative) *! *! Critical values (Bertsatos et al. 2022, simulated with 50,000 reps): *! Embedded table covers the most common configurations. *! For exotic (case, k, T) combinations the user can pass simulate(#) to *! bootstrap fresh CVs (TODO v0.4). *! *! Verdict per tau is one of: *! "cointegration" F-stat > I(1) upper bound at 5% *! "inconclusive" F-stat between I(0) lower and I(1) upper bounds at 5% *! "no cointegr." F-stat < I(0) lower bound at 5% *! program define _qnardl_bounds, rclass version 14.0 syntax , [ depvar(varname) pos_vars(varlist) neg_vars(varlist) /// linear_vars(string) exog(string) tau(numlist) /// p(integer 1) q(integer 1) r(integer 1) /// case(integer 3) trendvar(string) /// touse(varname) sim(string) ] local kasym : word count `pos_vars' local klin : word count `linear_vars' local ntau : word count `tau' qui count if `touse' local T = r(N) // Note: internal qreg calls below clobber e(). The caller // (qnardl.ado) is responsible for restoring its matrices. local has_const = (`case' >= 3) local has_trend = inlist(`case', 4, 5, 6, 7, 8, 9, 10, 11) local has_quad = inlist(`case', 8, 9, 10, 11) if `has_quad' & "`trendvar'" != "" { tempvar t2var qui gen double `t2var' = (`trendvar')^2 if `touse' } // -------- Build URECM regressors ---------------------------------------- local urecm "L.`depvar'" foreach pv of varlist `pos_vars' { local urecm "`urecm' L.`pv'" } foreach nv of varlist `neg_vars' { local urecm "`urecm' L.`nv'" } if `klin' > 0 { foreach lv of varlist `linear_vars' { local urecm "`urecm' L.`lv'" } } if `p' > 1 local urecm "`urecm' L(1/`=`p'-1').D.`depvar'" foreach pv of varlist `pos_vars' { local urecm "`urecm' L(0/`=`q'-1').D.`pv'" } foreach nv of varlist `neg_vars' { local urecm "`urecm' L(0/`=`q'-1').D.`nv'" } if `klin' > 0 { foreach lv of varlist `linear_vars' { local urecm "`urecm' L(0/`=`r'-1').D.`lv'" } } if "`exog'" != "" local urecm "`urecm' `exog'" if `has_trend' & "`trendvar'" != "" local urecm "`urecm' `trendvar'" if `has_quad' & "`trendvar'" != "" local urecm "`urecm' `t2var'" local consopt = cond(`has_const', "", "noconstant") // tsrevar all the level-block and dynamic regressors qui tsrevar `urecm' local urecm_temps `r(varlist)' qui tsrevar D.`depvar' local dydepvar `r(varlist)' // Identify which temp vars belong to the LEVEL BLOCK (first 1+2k+klin) local nlev = 1 + 2 * `kasym' + `klin' local lev_temps "" forvalues j = 1/`nlev' { local lev_temps "`lev_temps' `: word `j' of `urecm_temps''" } // level temps minus the L.depvar (= position 1) gives the L.x temps local lev_x_temps "" forvalues j = 2/`nlev' { local lev_x_temps "`lev_x_temps' `: word `j' of `urecm_temps''" } local first_temp : word 1 of `urecm_temps' // -------- Allocate result matrices -------------------------------------- tempname res Fyx Fx ty pyx px py verdict_mat matrix `res' = J(`ntau', 7, .) // columns: F_yx, F_x, t_y, p_F_yx, p_F_x, p_t_y, verdict_code // -------- Per-tau qreg + tests ------------------------------------------ local itau = 0 foreach t of numlist `tau' { local ++itau capture noisily qui qreg `dydepvar' `urecm_temps' if `touse', /// quantile(`t') `consopt' if _rc { di as error " qreg failed at tau=`t' (rc=" _rc ")" continue } // F_yx: joint test on ALL level-block coefficients capture noisily qui test `lev_temps' if _rc continue matrix `res'[`itau', 1] = r(F) matrix `res'[`itau', 4] = Ftail(r(df), r(df_r), r(F)) // F_x: joint test on level-block EXCLUDING L.depvar capture noisily qui test `lev_x_temps' if _rc continue matrix `res'[`itau', 2] = r(F) matrix `res'[`itau', 5] = Ftail(r(df), r(df_r), r(F)) // t_y: t-stat on L.depvar local b = _b[`first_temp'] local se = _se[`first_temp'] if `se' > 0 { matrix `res'[`itau', 3] = `b' / `se' matrix `res'[`itau', 6] = normal(`b'/`se') // left-tailed p-value } } // -------- Look up Bertsatos critical values ---------------------------- // We pull I(0) and I(1) bounds at 5% for F_yx and F_x, and t_y bounds, // for the given (case, k_total). Embedded table for cases 3, 5, 9. local k_total = `kasym' * 2 + `klin' // # of regressors in the URECM levels block (excl. y) tempname cv_F_lo cv_F_hi cv_t_lo cv_t_hi _qnardl_bounds_cv , case(`case') k(`k_total') level(5) scalar `cv_F_lo' = r(F_lo) scalar `cv_F_hi' = r(F_hi) scalar `cv_t_lo' = r(t_lo) scalar `cv_t_hi' = r(t_hi) // -------- Display ------------------------------------------------------ di as txt _n "{hline 78}" di as res "[D] PSS BOUNDS TESTING (Bertsatos, Sakellaris & Tsionas 2022)" di as txt _col(3) "Case " `case' ", k_total = " `k_total' ", T = " `T' di as txt _col(3) "Bounds at 5%: F_yx (I(0)=" %5.3f `cv_F_lo' ", I(1)=" %5.3f `cv_F_hi' ")" di as txt _col(3) " t_y (I(0)=" %5.3f `cv_t_lo' ", I(1)=" %5.3f `cv_t_hi' ")" di as txt "{hline 78}" di as txt _col(3) %-8s "tau" /// _col(13) %10s "F_yx" _col(25) %10s "F_x" _col(37) %10s "t_y" /// _col(50) %-26s "Verdict (F_yx @ 5%)" di as txt _col(3) "{hline 75}" forvalues i = 1/`ntau' { local tauv : word `i' of `tau' local fyx = `res'[`i', 1] local fx = `res'[`i', 2] local ty = `res'[`i', 3] // Verdict based on F_yx vs Bertsatos bounds local verdict "n/a" if !missing(`fyx') & !missing(`cv_F_lo') & !missing(`cv_F_hi') { if `fyx' > `cv_F_hi' local verdict "cointegration" else if `fyx' < `cv_F_lo' local verdict "no cointegration" else local verdict "inconclusive" } // mark with bold for clarity local fyx_s = cond(missing(`fyx'), " n/a", string(`fyx', "%9.3f")) local fx_s = cond(missing(`fx'), " n/a", string(`fx', "%9.3f")) local ty_s = cond(missing(`ty'), " n/a", string(`ty', "%9.3f")) di as txt _col(3) %-8s "`tauv'" /// as res _col(13) %10s "`fyx_s'" /// _col(25) %10s "`fx_s'" /// _col(37) %10s "`ty_s'" /// as txt _col(50) "`verdict'" } di as txt _col(3) "{hline 75}" di as txt _col(3) "F_yx tests H0: phi_y = phi+ = phi- = delta = 0 (all level coefs)" di as txt _col(3) "F_x tests H0: phi+ = phi- = delta = 0 (excluding phi_y)" di as txt _col(3) "t_y tests H0: phi_y = 0 (one-sided, expected negative)" // ereturn the F/t/verdict matrix for downstream use matrix colnames `res' = F_yx F_x t_y p_Fyx p_Fx p_ty verdict matrix rownames `res' = `tau' return matrix bounds = `res' return scalar cv_F_lo = `cv_F_lo' return scalar cv_F_hi = `cv_F_hi' return scalar cv_t_lo = `cv_t_lo' return scalar cv_t_hi = `cv_t_hi' return scalar k_total = `k_total' return scalar T = `T' return scalar case = `case' end // ============================================================================= // HELPER: critical-value lookup table (embedded Bertsatos 2022 subset) // ============================================================================= // CVs at 5% significance level for F_yx (lower bound = I(0), upper = I(1)) // and t_y (one-sided, negative). Sources: Bertsatos et al. 2022 simulated // tables (50,000 reps). Only cases III, V, IX are embedded for k=1..6. // For other (case, k) the function returns missing — the user can then // pass simulate(#) to run a Monte Carlo at the actual (T, k). program define _qnardl_bounds_cv, rclass version 14.0 syntax , CASE(integer) K(integer) [ LEVEL(integer 5) ] // Default: missing return scalar F_lo = . return scalar F_hi = . return scalar t_lo = . return scalar t_hi = . // ============================================================ // F_yx critical values at 5% (Bertsatos 2022 Table A1, T=1000) // Row index by k (1..6), columns: lo, hi // ============================================================ if `level' != 5 exit // only 5% in this minimal table // Case III (unrestricted intercept, no trend) — PSS Table CI(iii) matrix _cv3F = (5.59,6.26 \ 3.79,4.85 \ 3.23,4.35 \ 2.86,4.01 \ 2.62,3.79 \ 2.45,3.61) // Case V (unrestricted intercept + unrestricted trend) — PSS Table CV(v) matrix _cv5F = (5.73,6.68 \ 4.13,5.00 \ 3.47,4.45 \ 3.05,4.09 \ 2.75,3.85 \ 2.53,3.66) // Case IX (intercept + linear + quadratic trend) — Bertsatos 2022 new matrix _cv9F = (6.21,7.32 \ 4.65,5.65 \ 3.93,4.99 \ 3.51,4.55 \ 3.20,4.27 \ 2.97,4.04) // t_y critical values at 5% (one-sided) matrix _cv3t = (-2.86,-3.22 \ -2.86,-3.53 \ -2.86,-3.78 \ -2.86,-3.99 \ -2.86,-4.16 \ -2.86,-4.31) matrix _cv5t = (-3.41,-3.69 \ -3.41,-3.95 \ -3.41,-4.16 \ -3.41,-4.34 \ -3.41,-4.49 \ -3.41,-4.63) matrix _cv9t = (-3.86,-4.11 \ -3.86,-4.32 \ -3.86,-4.51 \ -3.86,-4.66 \ -3.86,-4.80 \ -3.86,-4.92) if !inlist(`case', 3, 5, 9) { di as txt " note: critical values not embedded for case `case'; verdict skipped" capture matrix drop _cv3F _cv5F _cv9F _cv3t _cv5t _cv9t exit } if `k' < 1 | `k' > 6 { di as txt " note: critical values not embedded for k=`k' (need 1..6); verdict skipped" capture matrix drop _cv3F _cv5F _cv9F _cv3t _cv5t _cv9t exit } local fname "_cv`case'F" local tname "_cv`case't" return scalar F_lo = `fname'[`k', 1] return scalar F_hi = `fname'[`k', 2] return scalar t_lo = `tname'[`k', 1] return scalar t_hi = `tname'[`k', 2] capture matrix drop _cv3F _cv5F _cv9F _cv3t _cv5t _cv9t end