*! surardl_estat 1.0.1 10oct2026 *! Postestimation for surardl *! Dr Merwan Roudane -- merwanroudane920@gmail.com *! github.com/merwanroudane program define surardl_estat, rclass version 14.0 if ("`e(cmd)'" != "surardl") { di as error "last estimates not found; surardl postestimation requires surardl" exit 301 } gettoken sub rest : 0, parse(" ,") local sub = lower("`sub'") if ("`sub'" == "bounds") _sa_estat_bounds `rest' else if ("`sub'" == "longrun") _sa_estat_longrun `rest' else if ("`sub'" == "shortrun") _sa_estat_shortrun `rest' else if ("`sub'" == "ecequal") _sa_estat_ecequal `rest' else if ("`sub'" == "lrequal") _sa_estat_lrequal `rest' else if ("`sub'" == "poolability") _sa_estat_pool `rest' else if ("`sub'" == "sigma") _sa_estat_sigma `rest' else if ("`sub'" == "serial") _sa_estat_serial `rest' else if ("`sub'" == "efficiency") _sa_estat_eff `rest' else if ("`sub'" == "graph") _sa_estat_graph `rest' else if ("`sub'" == "summary") _sa_estat_summary `rest' else { di as error "invalid estat subcommand `sub'" di as text "{p 4 4 2}Valid subcommands: {bf:bounds}, {bf:longrun}, {bf:shortrun}, {bf:ecequal}, {bf:lrequal}, {bf:poolability}, {bf:sigma}, {bf:serial}, {bf:efficiency}, {bf:graph}, {bf:summary}.{p_end}" exit 321 } return add end * =================================================================== * estat bounds * =================================================================== program define _sa_estat_bounds, rclass syntax [, BOOTstrap(integer 0) SEED(string) BOOTTarget(string) /// BOOTDgp(string) BOOTLags(integer 1) /// SIMulate(integer 0) SIMRho(real 0) noRULE ] if (`bootstrap' > 0 | `simulate' > 0) { * Refit the stored model with the requested resampling. The * original command line is reused so that any if/in and all * modelling choices are preserved, but every resampling option * is stripped first, otherwise they would be given twice. local cl `"`e(cmdline)'"' local cl = subinstr(`"`cl'"', "surardl ", "", 1) foreach drop in bootstrap boottarget bootdgp bootlags simulate /// simrho seed { _sa_stripopt `"`cl'"' `drop' local cl `"`r(cl)'"' } local cl = subinstr(`"`cl'"', "nodrift", "", .) if (strpos(`"`cl'"', ",") == 0) local cl `"`cl' ,"' if ("`seed'" == "") local seed = cond(e(seed) < ., e(seed), 20261010) if ("`boottarget'" == "") local boottarget equation if ("`bootdgp'" == "") local bootdgp rw di as text "{p 0 2 2}Refitting with the requested resampling; this replaces the stored results.{p_end}" local opt "" if (`bootstrap' > 0) local opt "`opt' bootstrap(`bootstrap') boottarget(`boottarget') bootdgp(`bootdgp') bootlags(`bootlags')" if (`simulate' > 0) local opt "`opt' simulate(`simulate') simrho(`simrho')" surardl `cl' seed(`seed') `opt' exit } _sa_display, noheader noctable nolrtable nodiagnostics tempname BB matrix `BB' = e(bounds) return matrix bounds = `BB' return scalar F_joint = e(F_joint) return scalar p_joint = e(p_joint) end * =================================================================== * estat longrun * =================================================================== program define _sa_estat_longrun, rclass syntax [, Level(cilevel) FIEller BOOTstrap(integer 0) SEED(string) ] local N = e(N_eq) local K = e(k_x) local xv "`e(indepvars)'" local eqn "`e(eqnames)'" tempname LR FI BC matrix `LR' = e(longrun) if (`level' != e(level)) { * recompute the normal interval at the requested level local z = invnormal(1 - (1 - `level'/100)/2) forvalues r = 1/`=`N'*`K'' { matrix `LR'[`r', 7] = `LR'[`r',3] - `z' * `LR'[`r',4] matrix `LR'[`r', 8] = `LR'[`r',3] + `z' * `LR'[`r',4] } } if ("`fieller'" != "") { mata: _sa_e_fieller(`level', "`FI'") } if (`bootstrap' > 0) { if ("`seed'" == "") local seed = cond(e(seed) < ., e(seed), 20261010) di as text "{p 0 2 2}Bootstrapping the long-run coefficients under the alternative ({bf:`bootstrap'} replications, seed {bf:`seed'})...{p_end}" mata: _sa_e_lrboot(`bootstrap', `seed', `level', "`BC'") } di "" di as text "{hline 80}" di as text "Long-run coefficients theta = -beta/lambda (eq. 82-84)" di as text "{hline 80}" di as text " Coef. Std. err. z P>|z| [" %2.0f `level' "% conf. int.]" di as text "{hline 80}" forvalues i = 1/`N' { local nm : word `i' of `eqn' di as text %-24s abbrev("`nm'", 24) forvalues k = 1/`K' { local r = (`i'-1)*`K' + `k' local x : word `k' of `xv' _sa_stars `=`LR'[`r',6]' di as text " " %-20s abbrev("`x'", 20) as result /// %11.4f `LR'[`r',3] " " %10.4f `LR'[`r',4] " " %7.2f `LR'[`r',5] /// " " %6.3f `LR'[`r',6] " " %10.4f `LR'[`r',7] " " %10.4f `LR'[`r',8] /// as text " `r(stars)'" if ("`fieller'" != "") { local ty = `FI'[`r',3] if (`ty' == 0) { di as text %-22s " Fieller set" as result /// _col(48) %10.4f `FI'[`r',1] " " %10.4f `FI'[`r',2] } else if (`ty' == 1) { di as text %-22s " Fieller set" as result _col(40) /// "unbounded: (-inf, " %6.3f `FI'[`r',1] "] U [" %6.3f `FI'[`r',2] ", +inf)" } else if (`ty' == 2) { di as text %-22s " Fieller set" as result _col(48) "entire real line" } else { di as text %-22s " Fieller set" as result _col(48) "empty" } } if (`bootstrap' > 0) { di as text %-22s " Bootstrap pct" as result /// _col(48) %10.4f `BC'[`r',2] " " %10.4f `BC'[`r',3] di as text %-22s " Bootstrap BC" as result /// _col(48) %10.4f `BC'[`r',4] " " %10.4f `BC'[`r',5] } } } di as text "{hline 80}" di as text "{p 0 2 2}Significance: * p<0.10, ** p<0.05, *** p<0.01. Standard errors by the delta method, eq. (85)-(91).{p_end}" if ("`fieller'" != "") { di as text "{p 0 2 2}The Fieller (1954) set inverts H0: beta + theta*lambda = 0 and never divides by lambda-hat. When lambda is weakly identified the set is unbounded, which is the honest answer; the delta-method interval cannot report that.{p_end}" } return matrix longrun = `LR' if ("`fieller'" != "") return matrix fieller = `FI' if (`bootstrap' > 0) return matrix bootci = `BC' end * =================================================================== * estat shortrun * =================================================================== program define _sa_estat_shortrun, rclass syntax [, Level(cilevel) ] local N = e(N_eq) local dep "`e(depvar)'" local xv "`e(indepvars)'" local eqn "`e(eqnames)'" tempname b V FIT matrix `b' = e(b) matrix `V' = e(V) matrix `FIT' = e(eqfit) local cn : colnames `b' local ce : coleq `b' di "" di as text "{hline 80}" di as text "Short-run (impact) coefficients and deterministic terms" di as text "{hline 80}" di as text " {c |} Coef. Std. err. z P>|z| [" %2.0f `level' "% conf. interval]" forvalues i = 1/`N' { local nm : word `i' of `eqn' local dfi = `FIT'[`i', 4] di as text "{hline 13}{c +}{hline 66}" di as text %-12s abbrev("`nm'", 12) " {c |}" local j 1 foreach c of local cn { local eq : word `j' of `ce' if ("`eq'" == "`nm'") { local islev 0 if ("`c'" == "L.`dep'") local islev 1 foreach x of local xv { if ("`c'" == "L.`x'") local islev 1 } if (`islev' == 0) { _sa_coefrow `b' `V' "`ce'" "`cn'" "`nm'" "`c'" `level' `dfi' "`e(small)'" } } local j = `j' + 1 } } di as text "{hline 13}{c BT}{hline 66}" di as text "{p 0 2 2}D.x is the impact elasticity, eq. (13c); _trend and _cons are the deterministic terms of the conditional ECM.{p_end}" end * =================================================================== * estat ecequal -- equality of the speeds of adjustment * =================================================================== program define _sa_estat_ecequal, rclass syntax [, MTest(string) ] _sa_padjcode "`mtest'" local pm = r(code) local pn "`r(name)'" local N = e(N_eq) local eqn "`e(eqnames)'" local dep "`e(depvar)'" tempname OUT JT mata: _sa_e_lamtest(`pm', "`OUT'", "`JT'") di "" di as text "{hline 80}" di as text "Equality of the speeds of adjustment H0: lambda_i = lambda_j (section 4.5c)" di as text "{hline 80}" di as text " Unit i Unit j chi2(1) p p(`pn')" di as text "{hline 80}" local r 1 forvalues i = 1/`=`N'-1' { forvalues j = `=`i'+1'/`N' { local ni : word `i' of `eqn' local nj : word `j' of `eqn' _sa_stars `=`OUT'[`r',5]' di as text " " %-22s abbrev("`ni'",22) " " %-22s abbrev("`nj'",22) /// as result %9.3f `OUT'[`r',3] " " %8.4f `OUT'[`r',4] " " /// %8.4f `OUT'[`r',5] as text " `r(stars)'" local r = `r' + 1 } } di as text "{hline 80}" di as text "Joint H0: lambda_1 = ... = lambda_`N' chi2(" %2.0f `JT'[1,2] ") = " /// as result %9.3f `JT'[1,1] as text " p = " as result %6.4f `JT'[1,3] di as text "{hline 80}" di as text "{p 0 2 2}Adjusted p-values use the `pn' correction over the " /// `=`N'*(`N'-1)/2' " pairwise tests. Section 4.5 of the paper reports a maximum pairwise rejection rate of 53% at N = 5, T = 50 without such a correction.{p_end}" di as text "{p 0 2 2}Ordering the units by lambda-hat reveals the adjustment-speed clusters discussed in section 5; see {bf:estat graph, lambda}.{p_end}" return matrix pairs = `OUT' return scalar chi2 = `JT'[1,1] return scalar df = `JT'[1,2] return scalar p = `JT'[1,3] end * =================================================================== * estat lrequal -- equality of long-run coefficients * =================================================================== program define _sa_estat_lrequal, rclass syntax [anything] [, MTest(string) ] local N = e(N_eq) local K = e(k_x) local xv "`e(indepvars)'" local eqn "`e(eqnames)'" _sa_padjcode "`mtest'" local pm = r(code) local pn "`r(name)'" local klist "" if ("`anything'" == "") { forvalues k = 1/`K' { local klist "`klist' `k'" } } else { foreach v of local anything { local pos 0 local j 1 foreach x of local xv { if ("`x'" == "`v'") local pos `j' local j = `j' + 1 } if (`pos' == 0) { di as error "`v' is not one of the long-run regressors (`xv')" exit 198 } local klist "`klist' `pos'" } } foreach k of local klist { local x : word `k' of `xv' tempname OUT JT mata: _sa_e_lrtest(`k', `pm', "`OUT'", "`JT'") di "" di as text "{hline 80}" di as text "Equality of long-run coefficients on `x' H0: theta_i = theta_j (Remark 10)" di as text "{hline 80}" di as text " Unit i Unit j difference chi2(1) p p(`pn')" di as text "{hline 80}" local r 1 forvalues i = 1/`=`N'-1' { forvalues j = `=`i'+1'/`N' { local ni : word `i' of `eqn' local nj : word `j' of `eqn' _sa_stars `=`OUT'[`r',6]' di as text " " %-22s abbrev("`ni'",22) " " %-18s abbrev("`nj'",18) /// as result %11.4f `OUT'[`r',3] " " %8.3f `OUT'[`r',4] " " /// %8.4f `OUT'[`r',5] " " %8.4f `OUT'[`r',6] as text " `r(stars)'" local r = `r' + 1 } } di as text "{hline 80}" di as text "Joint H0: theta_1 = ... = theta_`N' chi2(" %2.0f `JT'[1,2] ") = " /// as result %9.3f `JT'[1,1] as text " p = " as result %6.4f `JT'[1,3] di as text "{hline 80}" return matrix pairs_`k' = `OUT' } di as text "{p 0 2 2}Tests use the multivariate delta method of equations (92)-(93) with the full SUR covariance, so the cross-equation covariances that equation-by-equation OLS sets to zero are properly accounted for.{p_end}" end * =================================================================== * estat poolability * =================================================================== program define _sa_estat_pool, rclass syntax [, LAMbda SHORTrun ALL ] local N = e(N_eq) local K = e(k_x) local what 1 local lbl "long-run level coefficients beta" if ("`lambda'" != "") { local what 2 local lbl "speeds of adjustment lambda" } if ("`shortrun'" != "") { local what 3 local lbl "short-run coefficients" } if ("`all'" != "") { local what 4 local lbl "lambda and beta jointly" } tempname OUT mata: _sa_e_pool(`what', "`OUT'") di "" di as text "{hline 80}" di as text "Poolability test across the `N' equations (section 4.5a)" di as text "{hline 80}" di as text "H0: the `lbl' are equal across all units" if (`OUT'[1,1] >= .) { di as error "the restriction is not defined (unit-specific lag orders differ)" exit 498 } di as text " chi2(" %3.0f `OUT'[1,2] ") = " as result %10.4f `OUT'[1,1] /// as text " p = " as result %6.4f `OUT'[1,3] di as text "{hline 80}" di as text "{p 0 2 2}Simulation 4 of the paper finds this test over-rejects when N is large and T small (21% at N = 5, T = 50, nominal 5%); interpret cautiously below T = 200 for N = 5, or below T = 100 for N <= 3.{p_end}" return scalar chi2 = `OUT'[1,1] return scalar df = `OUT'[1,2] return scalar p = `OUT'[1,3] end * =================================================================== * estat sigma * =================================================================== program define _sa_estat_sigma, rclass syntax [, CORRelation ] local N = e(N_eq) local eqn "`e(eqnames)'" tempname S C SD BP matrix `S' = e(Sigma) matrix `C' = e(Corr) matrix `SD' = e(sigdiag) matrix `BP' = e(bpcd) di "" di as text "{hline 80}" di as text "Contemporaneous covariance matrix Sigma-hat (eq. 55-56, divisor: `e(sigmadf)')" di as text "{hline 80}" _sa_matprint `S', names(`eqn') format(%9.5f) di "" di as text "Cross-equation correlation matrix of the OLS residuals (eq. 97)" _sa_matprint `C', names(`eqn') format(%9.4f) di "" di as text "{hline 80}" di as text " Breusch-Pagan LM (eq. 96-99) chi2(" %3.0f `BP'[1,2] ") = " /// as result %10.4f `BP'[1,1] as text " p = " as result %6.4f `BP'[1,3] di as text " Pesaran CD (eq. 100) = " /// as result %10.4f `BP'[1,4] as text " p = " as result %6.4f `BP'[1,5] di as text " Mean absolute off-diagonal correlation = " as result %10.4f `BP'[1,6] di as text " T / N = " as result %10.3f `SD'[1,3] di as text " Reciprocal condition number = " as result %10.3e `SD'[1,4] di as text " Smallest eigenvalue = " as result %10.3e `SD'[1,5] di as text " Positive definite = " as result %10s cond(`SD'[1,6]==1,"yes","NO") di as text " Relative Frobenius norm (off-diagonal) = " as result %10.4f `SD'[1,8] di as text "{hline 80}" di as text "{p 0 2 2}By Proposition 1 the SUR estimator coincides with equation-by-equation OLS when Sigma is diagonal. The LM test is the primary diagnostic; the CD test is a robustness check with low power at small N.{p_end}" di as text "{p 0 2 2}The relative Frobenius norm is surardl's definition of the off-diagonal share of Sigma, ||Sigma - diag(Sigma)||_F / ||Sigma||_F. Simulation 7 of the paper reports a related quantity but does not define it.{p_end}" return matrix Sigma = `S' return matrix Corr = `C' return scalar bp = `BP'[1,1] return scalar bp_p = `BP'[1,3] return scalar cd = `BP'[1,4] return scalar cd_p = `BP'[1,5] end * =================================================================== * estat serial * =================================================================== program define _sa_estat_serial, rclass syntax [, LAGs(integer 1) ] if (`lags' < 1) { di as error "lags() must be at least 1" exit 198 } local N = e(N_eq) local eqn "`e(eqnames)'" tempname BG mata: st_matrix("`BG'", _surardl_bg(_surardl_reload(), st_matrix("e(resid)"), `lags')) di "" di as text "{hline 80}" di as text "Breusch-Godfrey LM test for serial correlation in the SUR residuals" di as text "H0: no autocorrelation up to order `lags' (Assumption 1(iii), Remark 5)" di as text "{hline 80}" di as text " LM df p" forvalues i = 1/`N' { local nm : word `i' of `eqn' local pb = `BG'[`i',3] local flag "" if (`pb' < 0.05) local flag " <- reject at 5%" di as text " " %-22s abbrev("`nm'",22) as result %10.4f `BG'[`i',1] /// " " %8.0f `BG'[`i',2] " " %10.4f `pb' as text "`flag'" } di as text "{hline 80}" di as text "{p 0 2 2}Rejection means the ARDL lag orders are too short. Raise {bf:lags()} in surardl until the conditional errors are white noise: only then is the contemporaneous Sigma, rather than a long-run covariance matrix, the correct GLS weight (Remark 4).{p_end}" return matrix bgodfrey = `BG' end * =================================================================== * estat efficiency -- eq. (95) * =================================================================== program define _sa_estat_eff, rclass local N = e(N_eq) local K = e(k_x) local dep "`e(depvar)'" local xv "`e(indepvars)'" local eqn "`e(eqnames)'" tempname EFF b matrix `EFF' = e(efficiency) matrix `b' = e(b) local cn : colnames `b' local ce : coleq `b' di "" di as text "{hline 80}" di as text "Efficiency gain from system estimation, 1 - SE(SUR)/SE(OLS) (eq. 95)" di as text "{hline 80}" di as text " SE(OLS) SE(SUR) gain" local j 1 local tot 0 local cnt 0 foreach c of local cn { local eq : word `j' of `ce' local show 0 if ("`c'" == "L.`dep'") local show 1 foreach x of local xv { if ("`c'" == "L.`x'") local show 1 } if (`show') { di as text " " %-14s abbrev("`eq'",14) " " %-10s abbrev("`c'",10) /// as result %10.5f `EFF'[`j',1] " " %11.5f `EFF'[`j',2] " " /// %9.1f 100*`EFF'[`j',3] as text "%" local tot = `tot' + `EFF'[`j',3] local cnt = `cnt' + 1 } local j = `j' + 1 } di as text "{hline 80}" di as text "Mean gain on the level terms: " as result %6.1f 100*`tot'/`cnt' as text "%" di as text "{p 0 2 2}Proposition 2 guarantees a non-negative asymptotic gain whenever Sigma is non-diagonal and the regressors differ across equations. Negative finite-sample entries are possible and simply reflect sampling error in Sigma-hat.{p_end}" return matrix efficiency = `EFF' end * =================================================================== * estat summary * =================================================================== program define _sa_estat_summary local N = e(N_eq) local K = e(k_x) local eqn "`e(eqnames)'" tempname FIT LAGS BND matrix `FIT' = e(eqfit) matrix `LAGS' = e(lags) matrix `BND' = e(bounds) di "" di as text "{hline 80}" di as text "Equation summary" di as text "{hline 80}" di as text " lags k df_r RSS(SUR) R2(OLS) F" forvalues i = 1/`N' { local nm : word `i' of `eqn' local lt = "`=`LAGS'[`i',1]'" forvalues k = 1/`K' { local lt "`lt',`=`LAGS'[`i',`k'+1]'" } di as text " " %-18s abbrev("`nm'",18) %-10s "`lt'" /// as result %7.0f `FIT'[`i',5] " " %7.0f `FIT'[`i',4] " " /// %11.5f `FIT'[`i',1] " " %10.4f `FIT'[`i',3] " " %9.3f `BND'[`i',1] } di as text "{hline 80}" end * =================================================================== * helper: square matrix printer with wrapped column blocks * =================================================================== program define _sa_matprint syntax anything(name=M) [, Names(string) Format(string) ] if ("`format'" == "") local format "%9.4f" local N = rowsof(`M') local nm "`names'" if ("`nm'" == "") { forvalues i = 1/`N' { local nm "`nm' r`i'" } } * 6 columns per block keeps every line inside 80 characters local per 6 local b0 1 while (`b0' <= `N') { local b1 = min(`b0' + `per' - 1, `N') di as text _col(16) _continue forvalues j = `b0'/`b1' { local cj : word `j' of `nm' di as text %10s abbrev("`cj'", 10) _continue } di "" di as text "{hline 15}{c +}" _continue di as text "{hline `=10*(`b1'-`b0'+1)'}" forvalues i = 1/`N' { local ci : word `i' of `nm' di as text %14s abbrev("`ci'", 14) " {c |}" _continue forvalues j = `b0'/`b1' { di as result `format' `M'[`i', `j'] " " _continue } di "" } local b0 = `b1' + 1 if (`b0' <= `N') di "" } end * =================================================================== * helper: remove one name(...) option from a command line, matching * the balanced closing parenthesis * =================================================================== program define _sa_stripopt, rclass args cl nm local out `"`cl'"' local p = strpos(lower(`"`out'"'), "`nm'(") while (`p' > 0) { local pre = substr(`"`out'"', 1, `p' - 1) local rest = substr(`"`out'"', `p' + length("`nm'") + 1, .) local d 1 local i 1 local n = length(`"`rest'"') while (`d' > 0 & `i' <= `n') { local ch = substr(`"`rest'"', `i', 1) if ("`ch'" == "(") local d = `d' + 1 if ("`ch'" == ")") local d = `d' - 1 local i = `i' + 1 } local out `"`pre'`=substr(`"`rest'"', `i', .)'"' local p = strpos(lower(`"`out'"'), "`nm'(") } return local cl `"`out'"' end * =================================================================== * helper: multiple-testing method code * =================================================================== program define _sa_padjcode, rclass args m local m = lower("`m'") if ("`m'" == "") local m holm if ("`m'" == "bonferroni") local c 1 else if ("`m'" == "holm") local c 2 else if ("`m'" == "sidak") local c 3 else if ("`m'" == "none") local c 0 else { di as error "mtest() must be bonferroni, holm, sidak or none" exit 198 } return scalar code = `c' return local name "`m'" end * =================================================================== * estat graph -- journal-style figures * =================================================================== program define _sa_estat_graph syntax [, LAMbda BOUnds LONGrun SIGma EFFiciency BOOTdist ALL /// Name(string) SAVing(string) VARiable(string) /// EQuation(integer 1) NODRAW SCHeme(string) ] local any "`lambda'`bounds'`longrun'`sigma'`efficiency'`bootdist'`all'" if ("`any'" == "") local all all if ("`scheme'" == "") local scheme "s2color" local gopt "graphregion(color(white) margin(medium)) plotregion(color(white)) scheme(`scheme')" local made "" if ("`lambda'" != "" | "`all'" != "") { _sa_g_lambda, `gopt' stub(`name') local made "`made' `r(nm)'" } if ("`bounds'" != "" | "`all'" != "") { _sa_g_bounds, `gopt' stub(`name') local made "`made' `r(nm)'" } if ("`longrun'" != "" | "`all'" != "") { _sa_g_longrun, `gopt' stub(`name') variable(`variable') local made "`made' `r(nm)'" } if ("`sigma'" != "" | "`all'" != "") { _sa_g_sigma, `gopt' stub(`name') local made "`made' `r(nm)'" } if ("`efficiency'" != "" | "`all'" != "") { _sa_g_eff, `gopt' stub(`name') local made "`made' `r(nm)'" } if ("`bootdist'" != "") { _sa_g_boot, `gopt' stub(`name') equation(`equation') local made "`made' `r(nm)'" } if ("`saving'" != "") { foreach g of local made { capture graph export "`saving'_`g'.png", name(`g') replace width(2000) } } end * ---- speed of adjustment, sorted: reveals the adjustment clusters ---- program define _sa_g_lambda, rclass syntax [, Stub(string) * ] local gname "`stub'lambda" local N = e(N_eq) local dep "`e(depvar)'" local eqn "`e(eqnames)'" tempname b V matrix `b' = e(b) matrix `V' = e(V) local z = invnormal(1 - (1 - e(level)/100)/2) preserve clear qui set obs `N' qui gen str32 eqname = "" qui gen double lam = . qui gen double lo = . qui gen double hi = . forvalues i = 1/`N' { local nmi : word `i' of `eqn' qui replace eqname = "`nmi'" in `i' local pos 0 local cn : colnames `b' local ce : coleq `b' local j 1 foreach c of local cn { local eq : word `j' of `ce' if ("`eq'" == "`nmi'" & "`c'" == "L.`dep'") local pos `j' local j = `j' + 1 } if (`pos' > 0) { local bb = `b'[1, `pos'] local ss = sqrt(`V'[`pos', `pos']) qui replace lam = `bb' in `i' qui replace lo = `bb' - `z'*`ss' in `i' qui replace hi = `bb' + `z'*`ss' in `i' } } sort lam qui gen int ord = _n qui levelsof ord, local(ol) local lbl "" forvalues i = 1/`N' { local nmi = eqname[`i'] local lbl `"`lbl' `i' "`=abbrev("`nmi'",14)'""' } twoway (rcap lo hi ord, horizontal lcolor(gs7) lwidth(medthin)) /// (scatter ord lam, msymbol(O) msize(medium) mcolor(navy)), /// yscale(reverse) ylabel(`lbl', angle(0) labsize(small) nogrid) /// xline(0, lcolor(cranberry) lpattern(dash) lwidth(thin)) /// ytitle("") xtitle("Speed of adjustment {&lambda}{subscript:i}") /// title("Error-correction speed of adjustment", size(medium)) /// subtitle("sorted; `=e(level)'% confidence intervals", size(small)) /// note("Stability requires {&lambda}{subscript:i} < 0 (section 2.5)." /// " Clusters of similar {&lambda} can be tested with estat ecequal.", size(vsmall)) /// legend(off) name(`gname', replace) `options' restore return local nm "`gname'" end * ---- bounds F statistics against the PSS critical bands ---- program define _sa_g_bounds, rclass syntax [, Stub(string) * ] local gname "`stub'bounds" local N = e(N_eq) local eqn "`e(eqnames)'" tempname BND CVF matrix `BND' = e(bounds) matrix `CVF' = e(cvF) local l10 = `CVF'[1,1] local u10 = `CVF'[1,2] local l05 = `CVF'[2,1] local u05 = `CVF'[2,2] preserve clear qui set obs `N' qui gen int ord = _n qui gen double F = . forvalues i = 1/`N' { qui replace F = `BND'[`i',1] in `i' } local lbl "" forvalues i = 1/`N' { local nmi : word `i' of `eqn' local lbl `"`lbl' `i' "`=abbrev("`nmi'",12)'""' } qui summ F local ymax = max(r(max)*1.15, `u05'*1.4) twoway (bar F ord, barwidth(0.55) color(navy%80)) /// (function y = `l05', range(0.4 `=`N'+0.6') lcolor(cranberry) /// lpattern(dash) lwidth(medthin)) /// (function y = `u05', range(0.4 `=`N'+0.6') lcolor(cranberry) /// lpattern(solid) lwidth(medthin)) /// (function y = `l10', range(0.4 `=`N'+0.6') lcolor(gs9) /// lpattern(shortdash) lwidth(vthin)) /// (function y = `u10', range(0.4 `=`N'+0.6') lcolor(gs9) /// lpattern(longdash) lwidth(vthin)), /// xlabel(`lbl', angle(45) labsize(small) nogrid) /// ylabel(, angle(0) labsize(small)) /// xtitle("") ytitle("SUR bounds F statistic") /// title("SUR-ARDL bounds test for a level relationship", size(medium)) /// subtitle("case `=e(case)', k = `=e(k_x)'; PSS (2001) bounds", size(small)) /// legend(order(2 "I(0) 5%" 3 "I(1) 5%" 4 "I(0) 10%" 5 "I(1) 10%") /// rows(1) size(vsmall) region(lstyle(none))) /// note("Above the I(1) line: reject H0. Below the I(0) line: do not reject. Between: inconclusive.", size(vsmall)) /// name(`gname', replace) `options' restore return local nm "`gname'" end * ---- forest plot of the long-run coefficients ---- program define _sa_g_longrun, rclass syntax [, Stub(string) VARiable(string) * ] local gname "`stub'longrun" local N = e(N_eq) local K = e(k_x) local xv "`e(indepvars)'" local eqn "`e(eqnames)'" local kk 1 if ("`variable'" != "") { local j 1 local kk 0 foreach x of local xv { if ("`x'" == "`variable'") local kk `j' local j = `j' + 1 } if (`kk' == 0) { di as error "`variable' is not one of the long-run regressors" exit 198 } } local xnm : word `kk' of `xv' tempname LR matrix `LR' = e(longrun) preserve clear qui set obs `N' qui gen int ord = _n qui gen double th = . qui gen double lo = . qui gen double hi = . forvalues i = 1/`N' { local r = (`i'-1)*`K' + `kk' qui replace th = `LR'[`r',3] in `i' qui replace lo = `LR'[`r',7] in `i' qui replace hi = `LR'[`r',8] in `i' } local lbl "" forvalues i = 1/`N' { local nmi : word `i' of `eqn' local lbl `"`lbl' `i' "`=abbrev("`nmi'",14)'""' } twoway (rcap lo hi ord, horizontal lcolor(gs7) lwidth(medthin)) /// (scatter ord th, msymbol(S) msize(medium) mcolor(dkgreen)), /// yscale(reverse) ylabel(`lbl', angle(0) labsize(small) nogrid) /// xline(0, lcolor(cranberry) lpattern(dash) lwidth(thin)) /// ytitle("") xtitle("Long-run coefficient on `xnm'") /// title("Long-run coefficients {&theta} = -{&beta}/{&lambda}", size(medium)) /// subtitle("delta-method `=e(level)'% intervals", size(small)) /// note("Delta-method intervals under-cover when |{&lambda}| is small; see estat longrun, fieller.", size(vsmall)) /// legend(off) name(`gname', replace) `options' restore return local nm "`gname'" end * ---- cross-equation correlation heatmap (base twoway only) ---- program define _sa_g_sigma, rclass syntax [, Stub(string) * ] local gname "`stub'sigma" local N = e(N_eq) local eqn "`e(eqnames)'" tempname C matrix `C' = e(Corr) preserve clear qui set obs `=`N'*`N'' qui gen int row = . qui gen int col = . qui gen double rho = . local z 1 forvalues i = 1/`N' { forvalues j = 1/`N' { qui replace row = `i' in `z' qui replace col = `j' in `z' qui replace rho = `C'[`i',`j'] in `z' local z = `z' + 1 } } qui gen str8 lab = string(rho, "%4.2f") local lbl "" forvalues i = 1/`N' { local nmi : word `i' of `eqn' local lbl `"`lbl' `i' "`=abbrev("`nmi'",10)'""' } * colour by sign and strength, using five discrete bands qui gen byte bnd = 1 if rho < -0.5 qui replace bnd = 2 if inrange(rho, -0.5, -0.2) qui replace bnd = 3 if inrange(rho, -0.2, 0.2) qui replace bnd = 4 if inrange(rho, 0.2, 0.5) qui replace bnd = 5 if rho > 0.5 local lay "" local cols "cranberry orange%50 gs14 eltblue navy" forvalues bb = 1/5 { local cc : word `bb' of `cols' local lay "`lay' (scatter row col if bnd==`bb', msymbol(S) msize(*4.6) mcolor(`cc') mlcolor(white) mlwidth(vthin))" } twoway `lay' (scatter row col, msymbol(none) mlabel(lab) /// mlabposition(0) mlabsize(vsmall) mlabcolor(black)), /// yscale(reverse) ylabel(`lbl', angle(0) labsize(small) nogrid) /// xlabel(`lbl', angle(45) labsize(small) nogrid) /// xscale(range(0.5 `=`N'+0.5')) yscale(range(0.5 `=`N'+0.5')) /// ytitle("") xtitle("") /// title("Cross-equation residual correlations", size(medium)) /// subtitle("Breusch-Pagan LM = `=string(e(bp),"%6.2f")', p = `=string(e(bp_p),"%5.3f")'", size(small)) /// legend(off) aspectratio(1) name(`gname', replace) `options' restore return local nm "`gname'" end * ---- efficiency gain on the level terms, eq. (95) ---- program define _sa_g_eff, rclass syntax [, Stub(string) * ] local gname "`stub'efficiency" local N = e(N_eq) local dep "`e(depvar)'" local xv "`e(indepvars)'" local eqn "`e(eqnames)'" tempname EFF b matrix `EFF' = e(efficiency) matrix `b' = e(b) local cn : colnames `b' local ce : coleq `b' preserve clear qui set obs `N' qui gen int ord = _n qui gen double g = . forvalues i = 1/`N' { local nmi : word `i' of `eqn' local j 1 foreach c of local cn { local eq : word `j' of `ce' if ("`eq'" == "`nmi'" & "`c'" == "L.`dep'") { qui replace g = 100*`EFF'[`j',3] in `i' } local j = `j' + 1 } } local lbl "" forvalues i = 1/`N' { local nmi : word `i' of `eqn' local lbl `"`lbl' `i' "`=abbrev("`nmi'",12)'""' } twoway (bar g ord, barwidth(0.55) color(dkgreen%80)), /// xlabel(`lbl', angle(45) labsize(small) nogrid) /// ylabel(, angle(0) labsize(small)) yline(0, lcolor(gs8)) /// xtitle("") ytitle("Reduction in SE of {&lambda} (%)") /// title("Efficiency gain from system estimation", size(medium)) /// subtitle("1 - SE(SUR)/SE(OLS), eq. (95)", size(small)) /// legend(off) name(`gname', replace) `options' restore return local nm "`gname'" end * ---- bootstrap null distribution with the observed statistic ---- program define _sa_g_boot, rclass syntax [, Stub(string) EQuation(integer 1) * ] local gname "`stub'bootdist" capture matrix list e(bootdist) if _rc { di as error "no bootstrap distribution stored; run surardl with bootstrap(#)" exit 498 } local N = e(N_eq) local eqn "`e(eqnames)'" if (!inrange(`equation', 1, `N')) { di as error "equation() must be between 1 and `N'" exit 198 } local nmi : word `equation' of `eqn' tempname BD BND BCV matrix `BD' = e(bootdist) matrix `BND' = e(bounds) matrix `BCV' = e(bootcv) local Fobs = `BND'[`equation',1] local c95 = `BCV'[`equation',2] preserve clear qui svmat double `BD', name(bs) qui keep bs`equation' qui rename bs`equation' Fstar qui drop if missing(Fstar) qui summ Fstar, detail local xmax = max(r(p99)*1.2, `Fobs'*1.15) twoway (histogram Fstar if Fstar <= `xmax', percent bin(40) /// color(navy%45) lcolor(white) lwidth(vvthin)) /// (pci 0 `Fobs' 12 `Fobs', lcolor(cranberry) lwidth(medium)) /// (pci 0 `c95' 12 `c95', lcolor(dkgreen) lpattern(dash) lwidth(medthin)), /// xscale(range(0 `xmax')) ylabel(, angle(0) labsize(small)) /// xtitle("Bootstrap F* under H0") ytitle("Percent") /// title("Bootstrap null distribution: `=abbrev("`nmi'",20)'", size(medium)) /// subtitle("`=e(reps)' replications, `=e(boottarget)' null, `=e(bootdgp)' marginal DGP", size(small)) /// legend(order(2 "observed F" 3 "bootstrap 5% CV") rows(1) /// size(vsmall) region(lstyle(none))) /// name(`gname', replace) `options' restore return local nm "`gname'" end * ------------------------------------------------------------------- * Display: journal-style tables built by hand (not ereturn display) * ------------------------------------------------------------------- program define _sa_display syntax [, Level(cilevel) noCTable noHEader noLRtable noBOUNDS /// noDIAGnostics ] local N = e(N_eq) local K = e(k_x) local T = e(T) local case = e(case) local dep "`e(depvar)'" local xv "`e(indepvars)'" local eqn "`e(eqnames)'" local small "`e(small)'" tempname LAGS FIT BND LR FI CVF CVT CVN DF DN DT SG BP BG BCV BPV S1 S0 EFF matrix `LAGS' = e(lags) matrix `FIT' = e(eqfit) matrix `BND' = e(bounds) matrix `LR' = e(longrun) matrix `FI' = e(fieller) matrix `CVF' = e(cvF) matrix `CVT' = e(cvt) matrix `CVN' = e(cvN) matrix `DF' = e(decF) matrix `DN' = e(decN) matrix `DT' = e(dect) matrix `SG' = e(sigdiag) matrix `BP' = e(bpcd) matrix `EFF' = e(efficiency) * =============================================================== * header * =============================================================== if ("`header'" == "") { di "" di as text "{hline 80}" di as text "Seemingly unrelated ARDL (SUR-ARDL) error-correction model" di as text "{hline 80}" local mname "Two-step feasible GLS (Zellner SUR)" if ("`e(method)'" == "iterated") local mname "Iterated FGLS (converges to ML)" if ("`e(method)'" == "ols") local mname "Equation-by-equation OLS (benchmark)" di as text %-42s "`mname'" _col(44) as text "Equations = " as result %8.0f `N' di as text %-42s "Dependent variable: D.`dep'" _col(44) as text "Time periods = " as result %8.0f `T' di as text %-42s "Case `case': `e(det)'" _col(44) as text "Observations = " as result %8.0f e(N) di as text %-42s "Panel: `e(ivar)' Time: `e(tvar)'" _col(44) as text "Parameters = " as result %8.0f e(M) local lagtxt = "ARDL(" + "`=`LAGS'[1,1]'" forvalues k = 1/`K' { local lagtxt "`lagtxt',`=`LAGS'[1,`k'+1]'" } local lagtxt "`lagtxt')" local homog 1 forvalues i = 2/`N' { forvalues j = 1/`=`K'+1' { if (`LAGS'[`i',`j'] != `LAGS'[1,`j']) local homog 0 } } if (`homog' == 0) local lagtxt "unit-specific (see e(lags))" if ("`e(lagsel)'" != "") local lagtxt "`lagtxt', selected by `=upper("`e(lagsel)'")'" di as text %-42s "Lag structure: `lagtxt'" _col(44) as text "Restrictions = " as result %8.0f e(df_q) if (e(n_gaps) > 0) { di as text "{p 0 2 2}Note: `=e(n_gaps)' complete period(s) outside the longest contiguous window were dropped.{p_end}" } if ("`e(method)'" == "iterated") { if (e(converged) == 0) { di as text "{p 0 2 2}Warning: iteration did not converge in `=e(iterations)' steps.{p_end}" } else { di as text "{p 0 2 2}Converged in `=e(iterations)' iterations.{p_end}" } } } * =============================================================== * coefficient table, grouped ADJ / LR / SR as in the EC form * =============================================================== if ("`ctable'" == "") { di "" di as text "{hline 80}" di as text "Conditional error-correction estimates " _col(44) "D.`dep' is the dependent variable" di as text "{hline 80}" di as text " {c |} Coef. Std. err. " _continue if ("`small'" != "") di as text "t P>|t| [" %2.0f `level' "% conf. interval]" else di as text "z P>|z| [" %2.0f `level' "% conf. interval]" local M = e(M) tempname b V matrix `b' = e(b) matrix `V' = e(V) local cn : colnames `b' local ce : coleq `b' forvalues i = 1/`N' { local nm : word `i' of `eqn' local dfi = `FIT'[`i', 4] di as text "{hline 13}{c +}{hline 66}" di as text %-12s abbrev("`nm'", 12) " {c |}" * ADJ block di as text " ADJ {c |}" _sa_coefrow `b' `V' "`ce'" "`cn'" "`nm'" "L.`dep'" `level' `dfi' "`small'" * LR block (lagged levels) di as text " LR {c |}" foreach x of local xv { _sa_coefrow `b' `V' "`ce'" "`cn'" "`nm'" "L.`x'" `level' `dfi' "`small'" } * SR block (everything else) di as text " SR {c |}" local j 1 foreach c of local cn { local eq : word `j' of `ce' if ("`eq'" == "`nm'") { local islev 0 if ("`c'" == "L.`dep'") local islev 1 foreach x of local xv { if ("`c'" == "L.`x'") local islev 1 } if (`islev' == 0) { _sa_coefrow `b' `V' "`ce'" "`cn'" "`nm'" "`c'" `level' `dfi' "`small'" } } local j = `j' + 1 } } di as text "{hline 13}{c BT}{hline 66}" di as text "{p 0 2 2}ADJ: speed of adjustment, eq. (8). LR: lagged level terms whose ratio to ADJ gives the long-run coefficients, eq. (82). SR: short-run dynamics and deterministic terms.{p_end}" } * =============================================================== * long-run coefficients * =============================================================== if ("`lrtable'" == "") { di "" di as text "{hline 80}" di as text "Long-run coefficients theta = -beta/lambda (delta method, eq. 82-91)" di as text "{hline 80}" di as text " Coef. Std. err. z P>|z| [" %2.0f `level' "% conf. int.]" di as text "{hline 80}" forvalues i = 1/`N' { local nm : word `i' of `eqn' di as text %-24s abbrev("`nm'", 24) forvalues k = 1/`K' { local r = (`i'-1)*`K' + `k' local x : word `k' of `xv' local th = `LR'[`r', 3] local se = `LR'[`r', 4] local zz = `LR'[`r', 5] local pp = `LR'[`r', 6] local lo = `LR'[`r', 7] local hi = `LR'[`r', 8] _sa_stars `pp' di as text " " %-20s abbrev("`x'", 20) as result /// %11.4f `th' " " %10.4f `se' " " %7.2f `zz' " " %6.3f `pp' /// " " %10.4f `lo' " " %10.4f `hi' as text " `r(stars)'" } } di as text "{hline 80}" di as text "{p 0 2 2}Significance: * p<0.10, ** p<0.05, *** p<0.01.{p_end}" di as text "{p 0 2 2}Simulation 5 of the paper documents severe under-coverage of the delta-method interval when |lambda| is small or T is moderate. Use {bf:estat longrun, fieller} for the Fieller (1954) set, which does not divide by lambda, or {bf:estat longrun, bootstrap(#)}.{p_end}" } * =============================================================== * bounds test * =============================================================== if ("`bounds'" == "") { local haveb = ("`e(boottarget)'" != "") if (`haveb') { matrix `BCV' = e(bootcv) matrix `BPV' = e(bootp) } local haves = (e(simreps) != . & e(simreps) > 0) if (`haves') { matrix `S1' = e(simcv1) matrix `S0' = e(simcv0) } di "" di as text "{hline 80}" di as text "SUR bounds test for a long-run level relationship H0: lambda_i = beta_i = 0" di as text "{hline 80}" di as text " F t " _continue if (`haveb') di as text "p(boot) decision (10%/5%/1%)" else di as text " decision (10%/5%/1%)" di as text "{hline 80}" forvalues i = 1/`N' { local nm : word `i' of `eqn' local Fi = `BND'[`i', 1] local ti = `BND'[`i', 3] local d10 = `DF'[`i', 1] local d05 = `DF'[`i', 2] local d01 = `DF'[`i', 4] _sa_dec `d10' local s10 "`r(d)'" _sa_dec `d05' local s05 "`r(d)'" _sa_dec `d01' local s01 "`r(d)'" if (`haveb') { di as text %-20s abbrev("`nm'", 20) as result %9.3f `Fi' " " %8.3f `ti' /// " " %7.3f `BPV'[`i',1] as text " `s10' `s05' `s01'" } else { di as text %-20s abbrev("`nm'", 20) as result %9.3f `Fi' " " %8.3f `ti' /// as text " `s10' `s05' `s01'" } } di as text "{hline 80}" di as text "Pesaran, Shin & Smith (2001) asymptotic bounds, case `case', k = `K':" di as text " 10% 5% 2.5% 1%" di as text " F " _continue forvalues j = 1/4 { di as result %7.2f `CVF'[`j',1] "/" %-6.2f `CVF'[`j',2] _continue } di "" di as text " t " _continue forvalues j = 1/4 { di as result %7.2f `CVT'[`j',1] "/" %-6.2f `CVT'[`j',2] _continue } di "" if (`CVN'[1,1] < .) { di as text "Narayan (2005) small-sample bounds, n = `T':" di as text " 10% 5% 1%" di as text " F " _continue forvalues j = 1/3 { di as result %7.2f `CVN'[`j',1] "/" %-6.2f `CVN'[`j',2] _continue } di "" } if (`haveb') { di as text "SUR residual-bootstrap critical values for F (`e(reps)' reps, seed `e(seed)', `e(boottarget)' null, `e(bootdgp)' marginal DGP):" di as text " 10% 5% 1% p-value" forvalues i = 1/`N' { local nm : word `i' of `eqn' di as text " " %-18s abbrev("`nm'", 18) as result /// %9.3f `BCV'[`i',1] " " %9.3f `BCV'[`i',2] " " %9.3f `BCV'[`i',3] /// " " %9.3f `BPV'[`i',1] } if (e(bfail) > 0) { di as text "{p 3 5 2}Note: `=e(bfail)' replication(s) failed and were discarded.{p_end}" } } if (`haves') { di as text "SUR simulated bounds for F (`e(simreps)' reps): I(0) lower / I(1) upper" di as text " 10% 5% 1%" forvalues i = 1/`N' { local nm : word `i' of `eqn' di as text " " %-18s abbrev("`nm'", 18) as result /// %7.2f `S0'[`i',1] "/" %-6.2f `S1'[`i',1] /// %7.2f `S0'[`i',2] "/" %-6.2f `S1'[`i',2] /// %7.2f `S0'[`i',3] "/" %-6.2f `S1'[`i',3] } } di as text "{hline 80}" di as text "Joint test across all `N' equations (eq. 78-79): " _continue di as result "F = " %7.3f e(F_joint) as text " chi2(" %2.0f e(df_joint) ") = " /// as result %8.3f e(chi2_joint) as text " p = " as result %6.4f e(p_joint) di as text "{p 0 2 2}Decision codes: {bf:R} reject (above the I(1) bound), {bf:-} inconclusive, {bf:N} no rejection (below the I(0) bound), {bf:.} bounds unavailable.{p_end}" if (`haveb' == 0 & `haves' == 0) { di as text "{p 0 2 2}{bf:Note:} the PSS and Narayan bounds above are calibrated for the {it:single-equation OLS} F statistic. The SUR statistic has a more dispersed null distribution in small samples, because Sigma is estimated jointly, so these bounds are liberal here. For SUR-specific values use {bf:bootstrap(}#{bf:)} (data-tailored, recommended) or {bf:simulate(}#{bf:)} (I(0)/I(1) bounds).{p_end}" } if (`N' >= 5 & `BP'[1,6] >= 0.5) { di as text "{p 0 2 2}{bf:Caution:} with N = `N' and mean |correlation| = " %4.2f `BP'[1,6] ", Simulation 2 of the paper puts the size distortion at its worst. Bootstrap critical values are not optional here.{p_end}" } if ("`e(fstat)'" != "wald") { di as text "{p 0 2 2}F reported from the restricted-system RSS form (eq. 73-76), fstat(`e(fstat)').{p_end}" } } * =============================================================== * diagnostics * =============================================================== if ("`diagnostics'" == "") { di "" di as text "{hline 80}" di as text "Diagnostics" di as text "{hline 80}" di as text "Cross-sectional dependence of the OLS residuals (justifies the SUR step):" di as text " Breusch-Pagan LM (eq. 96-99) chi2(" %3.0f `BP'[1,2] ") = " /// as result %9.4f `BP'[1,1] as text " p = " as result %6.4f `BP'[1,3] di as text " Pesaran CD (eq. 100) = " /// as result %9.4f `BP'[1,4] as text " p = " as result %6.4f `BP'[1,5] di as text " Mean absolute off-diagonal correlation = " as result %9.4f `BP'[1,6] if (`BP'[1,3] > 0.10) { di as text "{p 3 5 2}The LM test does not reject a diagonal Sigma. By Proposition 1 the SUR estimator then offers no efficiency gain over equation-by-equation OLS.{p_end}" } di as text "Covariance matrix Sigma-hat (`=`SG'[1,1]' x `=`SG'[1,1]'):" di as text " T/N ratio = " as result %9.3f `SG'[1,3] di as text " Reciprocal condition number = " as result %9.3e `SG'[1,4] di as text " Positive definite = " as result %9s cond(`SG'[1,6]==1, "yes", "NO") di as text " Relative Frobenius norm of off-diagonal = " as result %9.4f `SG'[1,8] if (`SG'[1,3] < 5) { di as text "{p 3 5 2}Warning: T/N < 5. Section 2.3 of the paper warns that Sigma-hat is estimated imprecisely in this regime.{p_end}" } capture matrix `BG' = e(bgodfrey) if (_rc == 0) { di as text "Breusch-Godfrey LM test on the SUR residuals (Assumption 1(iii), Remark 5):" di as text " LM df p" forvalues i = 1/`N' { local nm : word `i' of `eqn' local pb = `BG'[`i',3] local flag "" if (`pb' < 0.05) local flag " <- serial correlation at 5%" di as text " " %-18s abbrev("`nm'", 18) as result %10.4f `BG'[`i',1] /// " " %7.0f `BG'[`i',2] " " %8.4f `pb' as text "`flag'" } di as text "{p 3 5 2}If the null is rejected, raise the ARDL lag orders until the conditional errors are serially uncorrelated; Sigma rather than a long-run covariance is the correct weight only then (Remark 4).{p_end}" } * speed of adjustment sign check local nneg 0 forvalues i = 1/`N' { tempname bb matrix `bb' = e(b) local nm : word `i' of `eqn' local lam = _b[`nm':L.`dep'] if (`lam' >= 0) { local nneg = `nneg' + 1 di as text "{p 3 5 2}Warning: lambda = " %6.3f `lam' " >= 0 for `nm'; the error-correction mechanism is not stable, so no valid long-run relationship exists for that unit whatever the F statistic shows (section 2.5).{p_end}" } } di as text "{hline 80}" } end program define _sa_coefrow args b V ce cn nm cname level dfr small local j 1 local pos 0 foreach c of local cn { local eq : word `j' of `ce' if ("`eq'" == "`nm'" & "`c'" == "`cname'") local pos `j' local j = `j' + 1 } if (`pos' == 0) exit local bb = `b'[1, `pos'] local vv = `V'[`pos', `pos'] if (`vv' <= 0) exit local se = sqrt(`vv') local zz = `bb' / `se' if ("`small'" != "") { local pp = 2 * ttail(`dfr', abs(`zz')) local cv = invttail(`dfr', (1 - `level'/100)/2) } else { local pp = 2 * normal(-abs(`zz')) local cv = invnormal(1 - (1 - `level'/100)/2) } local lo = `bb' - `cv' * `se' local hi = `bb' + `cv' * `se' _sa_stars `pp' di as text %12s abbrev("`cname'", 12) " {c |}" as result /// %11.4f `bb' " " %11.4f `se' " " %7.2f `zz' " " %6.3f `pp' " " /// %10.4f `lo' " " %10.4f `hi' as text " `r(stars)'" end program define _sa_stars, rclass args p local s " " if (`p' < 0.10) local s "* " if (`p' < 0.05) local s "** " if (`p' < 0.01) local s "***" return local stars "`s'" end program define _sa_dec, rclass args d local s "." if (`d' == 0) local s "N" if (`d' == 1) local s "-" if (`d' == 2) local s "R" return local d "`s'" end * ------------------------------------------------------------------- * Mata library. * * Stata compiles this block into a namespace private to this ado-file, * so the same library is carried by surardl.ado and surardl_estat.ado. * Edit the _sa_*.mata sources and re-run build.py; do not edit here. * ------------------------------------------------------------------- version 14.0 mata: // ================================================================ // from _sa_cv.mata // ================================================================ // ------------------------------------------------------------------ // Pesaran, Shin & Smith (2001), Table CI: asymptotic critical value // bounds for the bounds F statistic. Rows are k = 0..10; the returned // 4 x 2 matrix has rows 10%, 5%, 2.5%, 1% and columns I(0), I(1). // ------------------------------------------------------------------ real matrix _surardl_cvF_pss(real scalar cs, real scalar k) { real matrix T real rowvector r if (cs == 1) { T = ( 3, 3, 4.2, 4.2, 5.47, 5.47, 7.17, 7.17 \ 2.44, 3.28, 3.15, 4.11, 3.88, 4.92, 4.81, 6.02 \ 2.17, 3.19, 2.72, 3.83, 3.22, 4.5, 3.88, 5.3 \ 2.01, 3.1, 2.45, 3.63, 2.87, 4.16, 3.42, 4.84 \ 1.9, 3.01, 2.26, 3.48, 2.62, 3.9, 3.07, 4.44 \ 1.81, 2.93, 2.14, 3.34, 2.44, 3.71, 2.82, 4.21 \ 1.75, 2.87, 2.04, 3.24, 2.32, 3.59, 2.66, 4.05 \ 1.7, 2.83, 1.97, 3.18, 2.22, 3.49, 2.54, 3.91 \ 1.66, 2.79, 1.91, 3.11, 2.15, 3.4, 2.45, 3.79 \ 1.63, 2.75, 1.86, 3.05, 2.08, 3.33, 2.34, 3.68 \ 1.6, 2.72, 1.82, 2.99, 2.02, 3.27, 2.26, 3.6 ) } if (cs == 2) { T = ( 3.8, 3.8, 4.6, 4.6, 5.39, 5.39, 6.44, 6.44 \ 3.02, 3.51, 3.62, 4.16, 4.18, 4.79, 4.94, 5.58 \ 2.63, 3.35, 3.1, 3.87, 3.55, 4.38, 4.13, 5 \ 2.37, 3.2, 2.79, 3.67, 3.15, 4.08, 3.65, 4.66 \ 2.2, 3.09, 2.56, 3.49, 2.88, 3.87, 3.29, 4.37 \ 2.08, 3, 2.39, 3.38, 2.7, 3.73, 3.06, 4.15 \ 1.99, 2.94, 2.27, 3.28, 2.55, 3.61, 2.88, 3.99 \ 1.92, 2.89, 2.17, 3.21, 2.43, 3.51, 2.73, 3.9 \ 1.85, 2.85, 2.11, 3.15, 2.33, 3.42, 2.62, 3.77 \ 1.8, 2.8, 2.04, 3.08, 2.24, 3.35, 2.5, 3.68 \ 1.76, 2.77, 1.98, 3.04, 2.18, 3.28, 2.41, 3.61 ) } if (cs == 3) { T = ( 6.58, 6.58, 8.21, 8.21, 9.8, 9.8, 11.79, 11.79 \ 4.04, 4.78, 4.94, 5.73, 5.77, 6.68, 6.84, 7.84 \ 3.17, 4.14, 3.79, 4.85, 4.41, 5.52, 5.15, 6.36 \ 2.72, 3.77, 3.23, 4.35, 3.69, 4.89, 4.29, 5.61 \ 2.45, 3.52, 2.86, 4.01, 3.25, 4.49, 3.74, 5.06 \ 2.26, 3.35, 2.62, 3.79, 2.96, 4.18, 3.41, 4.68 \ 2.12, 3.23, 2.45, 3.61, 2.75, 3.99, 3.15, 4.43 \ 2.03, 3.13, 2.32, 3.5, 2.6, 3.84, 2.96, 4.26 \ 1.95, 3.06, 2.22, 3.39, 2.48, 3.7, 2.79, 4.1 \ 1.88, 2.99, 2.14, 3.3, 2.37, 3.6, 2.65, 3.97 \ 1.83, 2.94, 2.06, 3.24, 2.28, 3.5, 2.54, 3.86 ) } if (cs == 4) { T = ( 5.37, 5.37, 6.29, 6.29, 7.14, 7.14, 8.26, 8.26 \ 4.05, 4.49, 4.68, 5.15, 5.3, 5.83, 6.1, 6.73 \ 3.38, 4.02, 3.88, 4.61, 4.37, 5.16, 4.99, 5.85 \ 2.97, 3.74, 3.38, 4.23, 3.8, 4.68, 4.3, 5.23 \ 2.68, 3.53, 3.05, 3.97, 3.4, 4.36, 3.81, 4.92 \ 2.49, 3.38, 2.81, 3.76, 3.11, 4.13, 3.5, 4.63 \ 2.33, 3.25, 2.63, 3.62, 2.9, 3.94, 3.27, 4.39 \ 2.22, 3.17, 2.5, 3.5, 2.76, 3.81, 3.07, 4.23 \ 2.13, 3.09, 2.38, 3.41, 2.62, 3.7, 2.93, 4.06 \ 2.05, 3.02, 2.3, 3.33, 2.52, 3.6, 2.79, 3.93 \ 1.98, 2.97, 2.21, 3.25, 2.42, 3.52, 2.68, 3.84 ) } if (cs == 5) { T = ( 9.81, 9.81, 11.64, 11.64, 13.36, 13.36, 15.73, 15.73 \ 5.59, 6.26, 6.56, 7.3, 7.46, 8.27, 8.74, 9.63 \ 4.19, 5.06, 4.87, 5.85, 5.49, 6.59, 6.34, 7.52 \ 3.47, 4.45, 4.01, 5.07, 4.52, 5.62, 5.17, 6.36 \ 3.03, 4.06, 3.47, 4.57, 3.89, 5.07, 4.4, 5.72 \ 2.75, 3.79, 3.12, 4.25, 3.47, 4.67, 3.93, 5.23 \ 2.53, 3.59, 2.87, 4, 3.19, 4.38, 3.6, 4.9 \ 2.38, 3.45, 2.69, 3.83, 2.98, 4.16, 3.34, 4.63 \ 2.26, 3.34, 2.55, 3.68, 2.82, 4.02, 3.15, 4.43 \ 2.16, 3.24, 2.43, 3.56, 2.67, 3.87, 2.97, 4.24 \ 2.07, 3.16, 2.33, 3.46, 2.56, 3.76, 2.84, 4.1 ) } if (k < 0 | k > 10) { return(J(4, 2, .)) } r = T[k + 1, .] return((r[1], r[2] \ r[3], r[4] \ r[5], r[6] \ r[7], r[8])) } // ------------------------------------------------------------------ // Pesaran, Shin & Smith (2001), Table CII: asymptotic critical value // bounds for the bounds t statistic on lambda. PSS tabulate cases I, // III and V only. Because the deterministic restrictions are ignored // for this single-coefficient test (Kripfganz & Schneider 2020, p. // 1461), case II uses the case III bounds and case IV the case V ones. // ------------------------------------------------------------------ real matrix _surardl_cvt_pss(real scalar cs, real scalar k) { real matrix T real rowvector r real scalar c c = cs if (cs == 2) { c = 3 } if (cs == 4) { c = 5 } if (c == 1) { T = ( -1.62, -1.62, -1.95, -1.95, -2.24, -2.24, -2.58, -2.58 \ -1.62, -2.28, -1.95, -2.6, -2.24, -2.9, -2.58, -3.22 \ -1.62, -2.68, -1.95, -3.02, -2.24, -3.31, -2.58, -3.66 \ -1.62, -3, -1.95, -3.33, -2.24, -3.64, -2.58, -3.97 \ -1.62, -3.26, -1.95, -3.6, -2.24, -3.89, -2.58, -4.23 \ -1.62, -3.49, -1.95, -3.83, -2.24, -4.12, -2.58, -4.44 \ -1.62, -3.7, -1.95, -4.04, -2.24, -4.34, -2.58, -4.67 \ -1.62, -3.9, -1.95, -4.23, -2.24, -4.54, -2.58, -4.88 \ -1.62, -4.09, -1.95, -4.43, -2.24, -4.72, -2.58, -5.07 \ -1.62, -4.26, -1.95, -4.61, -2.24, -4.89, -2.58, -5.25 \ -1.62, -4.42, -1.95, -4.76, -2.24, -5.06, -2.58, -5.44 ) } if (c == 3) { T = ( -2.57, -2.57, -2.86, -2.86, -3.13, -3.13, -3.43, -3.43 \ -2.57, -2.91, -2.86, -3.22, -3.13, -3.5, -3.43, -3.82 \ -2.57, -3.21, -2.86, -3.53, -3.13, -3.8, -3.43, -4.1 \ -2.57, -3.46, -2.86, -3.78, -3.13, -4.05, -3.43, -4.37 \ -2.57, -3.66, -2.86, -3.99, -3.13, -4.26, -3.43, -4.6 \ -2.57, -3.86, -2.86, -4.19, -3.13, -4.46, -3.43, -4.79 \ -2.57, -4.04, -2.86, -4.38, -3.13, -4.66, -3.43, -4.99 \ -2.57, -4.23, -2.86, -4.57, -3.13, -4.85, -3.43, -5.19 \ -2.57, -4.4, -2.86, -4.72, -3.13, -5.02, -3.43, -5.37 \ -2.57, -4.56, -2.86, -4.88, -3.13, -5.18, -3.42, -5.54 \ -2.57, -4.69, -2.86, -5.03, -3.13, -5.34, -3.43, -5.68 ) } if (c == 5) { T = ( -3.13, -3.13, -3.41, -3.41, -3.65, -3.66, -3.96, -3.97 \ -3.13, -3.4, -3.41, -3.69, -3.65, -3.96, -3.96, -4.26 \ -3.13, -3.63, -3.41, -3.95, -3.65, -4.2, -3.96, -4.53 \ -3.13, -3.84, -3.41, -4.16, -3.65, -4.42, -3.96, -4.73 \ -3.13, -4.04, -3.41, -4.36, -3.65, -4.62, -3.96, -4.96 \ -3.13, -4.21, -3.41, -4.52, -3.65, -4.79, -3.96, -5.13 \ -3.13, -4.37, -3.41, -4.69, -3.65, -4.96, -3.96, -5.31 \ -3.13, -4.53, -3.41, -4.85, -3.65, -5.14, -3.96, -5.49 \ -3.13, -4.68, -3.41, -5.01, -3.65, -5.3, -3.96, -5.65 \ -3.13, -4.82, -3.41, -5.15, -3.65, -5.44, -3.96, -5.79 \ -3.13, -4.96, -3.41, -5.29, -3.65, -5.59, -3.96, -5.94 ) } if (k < 0 | k > 10) { return(J(4, 2, .)) } r = T[k + 1, .] return((r[1], r[2] \ r[3], r[4] \ r[5], r[6] \ r[7], r[8])) } // ------------------------------------------------------------------ // Narayan (2005) small-sample critical value bounds for the bounds F // statistic: cases II-V, effective sample size n = 30(5)80, k = 0..7, // at the 10%, 5% and 1% levels, linearly interpolated in n. Column 1 // of each stored table is n; the remaining 48 columns run level-major // (10%, 5%, 1%), then k = 0..7, then (I(0), I(1)). // // Two documented departures from the printed tables: // (a) case II, 10%, n = 80: the published I(0) entry for k = 5 // repeats the k = 4 entry (2.303). It is replaced by 2.180, the // value printed for the same row on Narayan's case IV page. // (b) case IV, n = 80: all three published rows are verbatim copies // of the case II n = 80 rows and are evidently a typesetting // error. They are not stored, so case IV interpolation is // clamped at n = 75, which is mildly conservative. // ------------------------------------------------------------------ real matrix _surardl_cvF_nar(real scalar cs, real scalar k, real scalar n) { real matrix T, O real scalar i, nr, lo, hi, w, j real rowvector a, b if (cs < 2 | cs > 5 | k < 0 | k > 7) { return(J(3, 2, .)) } if (cs == 2) { T = J(11, 49, .) T[1, .] = (30, 4.025, 4.025, 3.303, 3.797, 2.915, 3.695, 2.676, 3.586, 2.525, 3.56, 2.407, 3.517, 2.334, 3.515, 2.277, 3.498, 5.07, 5.07, 4.09, 4.663, 3.538, 4.428, 3.272, 4.306, 3.058, 4.223, 2.91, 4.193, 2.794, 4.148, 2.73, 4.163, 7.595, 7.595, 6.027, 6.76, 5.155, 6.265, 4.614, 5.966, 4.28, 5.84, 4.134, 5.761, 3.976, 5.691, 3.864, 5.694) T[2, .] = (35, 3.98, 3.98, 3.223, 3.757, 2.845, 3.623, 2.618, 3.532, 2.46, 3.46, 2.331, 3.417, 2.254, 3.388, 2.196, 3.37, 4.945, 4.945, 3.957, 4.53, 3.478, 4.335, 3.164, 4.194, 2.947, 4.088, 2.804, 4.013, 2.685, 3.96, 2.597, 3.907, 7.35, 7.35, 5.763, 6.48, 4.948, 6.028, 4.428, 5.816, 4.093, 5.532, 3.9, 5.419, 3.713, 5.326, 3.599, 5.23) T[3, .] = (40, 3.955, 3.955, 3.21, 3.73, 2.835, 3.585, 2.592, 3.454, 2.427, 3.395, 2.306, 3.353, 2.218, 3.314, 2.152, 3.296, 4.96, 4.96, 3.937, 4.523, 3.435, 4.26, 3.1, 4.088, 2.893, 4, 2.734, 3.92, 2.618, 3.863, 2.523, 3.829, 7.22, 7.22, 5.593, 6.333, 4.77, 5.855, 4.31, 5.544, 3.967, 5.455, 3.657, 5.256, 3.505, 5.121, 3.402, 5.031) T[4, .] = (45, 3.95, 3.95, 3.19, 3.73, 2.788, 3.54, 2.56, 3.428, 2.402, 3.345, 2.276, 3.297, 2.188, 3.254, 2.131, 3.223, 4.895, 4.895, 3.877, 4.46, 3.368, 4.203, 3.078, 4.022, 2.85, 3.905, 2.694, 3.829, 2.591, 3.766, 2.504, 3.723, 7.265, 7.265, 5.607, 6.193, 4.8, 5.725, 4.27, 5.412, 3.892, 5.173, 3.674, 5.019, 3.54, 4.931, 3.383, 4.832) T[5, .] = (50, 3.935, 3.935, 3.177, 3.653, 2.788, 3.513, 2.538, 3.398, 2.372, 3.32, 2.259, 3.264, 2.17, 3.22, 2.099, 3.181, 4.815, 4.815, 3.86, 4.44, 3.368, 4.178, 3.048, 4.002, 2.823, 3.872, 2.67, 3.781, 2.55, 3.708, 2.457, 3.65, 7.065, 7.065, 5.503, 6.24, 4.695, 5.758, 4.188, 5.328, 3.845, 5.15, 3.593, 4.981, 3.424, 4.88, 3.282, 4.73) T[6, .] = (55, 3.9, 3.9, 3.143, 3.67, 2.748, 3.495, 2.508, 3.356, 2.345, 3.28, 2.226, 3.241, 2.139, 3.204, 2.069, 3.148, 4.795, 4.795, 3.79, 4.393, 3.303, 4.1, 2.982, 3.942, 2.763, 3.813, 2.617, 3.743, 2.49, 3.658, 2.414, 3.608, 6.965, 6.965, 5.377, 6.047, 4.61, 5.563, 4.118, 5.2, 3.738, 4.947, 3.543, 4.839, 3.33, 4.708, 3.194, 4.562) T[7, .] = (60, 3.88, 3.88, 3.127, 3.65, 2.738, 3.465, 2.496, 3.346, 2.323, 3.273, 2.204, 3.21, 2.114, 3.153, 2.044, 3.104, 4.78, 4.78, 3.803, 4.363, 3.288, 4.07, 2.962, 3.91, 2.743, 3.792, 2.589, 3.683, 2.456, 3.598, 2.373, 3.54, 6.96, 6.96, 5.383, 6.033, 4.558, 5.59, 4.068, 5.25, 3.71, 4.965, 3.451, 4.764, 3.293, 4.615, 3.129, 4.507) T[8, .] = (65, 3.88, 3.88, 3.143, 3.623, 2.74, 3.455, 2.492, 3.35, 2.335, 3.252, 2.209, 3.201, 2.12, 3.145, 2.043, 3.094, 4.78, 4.78, 3.787, 4.343, 3.285, 4.07, 2.976, 3.896, 2.75, 3.755, 2.596, 3.677, 2.473, 3.583, 2.373, 3.519, 6.825, 6.825, 5.35, 6.017, 4.538, 5.475, 4.056, 5.158, 3.725, 4.94, 3.43, 4.721, 3.225, 4.571, 3.092, 4.478) T[9, .] = (70, 3.875, 3.875, 3.12, 3.623, 2.73, 3.445, 2.482, 3.31, 2.32, 3.232, 2.193, 3.161, 2.1, 3.121, 2.024, 3.079, 4.75, 4.75, 3.78, 4.327, 3.243, 4.043, 2.924, 3.86, 2.725, 3.718, 2.564, 3.65, 2.451, 3.559, 2.351, 3.498, 6.74, 6.74, 5.157, 5.957, 4.398, 5.463, 3.916, 5.088, 3.608, 4.86, 3.373, 4.717, 3.18, 4.596, 3.034, 4.426) T[10, .] = (75, 3.895, 3.895, 3.133, 3.597, 2.725, 3.455, 2.482, 3.334, 2.313, 3.228, 2.196, 3.166, 2.103, 3.111, 2.023, 3.068, 4.76, 4.76, 3.777, 4.32, 3.253, 4.065, 2.946, 3.862, 2.725, 3.718, 2.574, 3.641, 2.449, 3.55, 2.36, 3.478, 6.915, 6.915, 5.26, 5.957, 4.458, 5.41, 4.048, 5.092, 3.687, 4.842, 3.427, 4.62, 3.219, 4.526, 3.057, 4.413) T[11, .] = (80, 3.87, 3.87, 3.113, 3.61, 2.713, 3.453, 2.474, 3.312, 2.303, 3.22, 2.18, 3.154, 2.088, 3.103, 2.017, 3.052, 4.725, 4.725, 3.74, 4.303, 3.235, 4.053, 2.92, 3.838, 2.688, 3.698, 2.55, 3.606, 2.431, 3.518, 2.336, 3.458, 6.695, 6.695, 5.157, 5.917, 4.358, 5.393, 3.908, 5.004, 3.602, 4.787, 3.351, 4.587, 3.173, 4.485, 3.021, 4.35) } if (cs == 3) { T = J(11, 49, .) T[1, .] = (30, 6.84, 6.84, 4.29, 5.08, 3.437, 4.47, 3.008, 4.15, 2.752, 3.994, 2.578, 3.858, 2.457, 3.797, 2.384, 3.728, 8.77, 8.77, 5.395, 6.35, 4.267, 5.473, 3.71, 5.018, 3.354, 4.774, 3.125, 4.608, 2.97, 4.499, 2.875, 4.445, 13.68, 13.68, 8.17, 9.285, 6.183, 7.873, 5.333, 7.063, 4.768, 6.67, 4.537, 6.37, 4.27, 6.211, 4.104, 6.151) T[2, .] = (35, 6.81, 6.81, 4.225, 5.05, 3.393, 4.41, 2.958, 4.1, 2.696, 3.898, 2.508, 3.763, 2.387, 3.671, 2.3, 3.606, 8.64, 8.64, 5.29, 6.175, 4.183, 5.333, 3.615, 4.913, 3.276, 4.63, 3.037, 4.443, 2.864, 4.324, 2.753, 4.209, 13.29, 13.29, 7.87, 8.96, 6.14, 7.607, 5.198, 6.845, 4.59, 6.368, 4.257, 6.04, 4.016, 5.797, 3.841, 5.686) T[3, .] = (40, 6.76, 6.76, 4.235, 5, 3.373, 4.377, 2.933, 4.02, 2.66, 3.838, 2.483, 3.708, 2.353, 3.599, 2.26, 3.534, 8.57, 8.57, 5.26, 6.16, 4.133, 5.26, 3.548, 4.803, 3.202, 4.544, 2.962, 4.338, 2.797, 4.211, 2.676, 4.13, 13.07, 13.07, 7.625, 8.825, 5.893, 7.337, 5.018, 6.61, 4.428, 6.25, 4.045, 5.898, 3.8, 5.643, 3.644, 5.464) T[4, .] = (45, 6.76, 6.76, 4.225, 5.02, 3.33, 4.347, 2.893, 3.983, 2.638, 3.772, 2.458, 3.647, 2.327, 3.541, 2.238, 3.461, 8.59, 8.59, 5.235, 6.135, 4.083, 5.207, 3.535, 4.733, 3.178, 4.45, 2.922, 4.268, 2.764, 4.123, 2.643, 4.004, 12.93, 12.93, 7.74, 8.65, 5.92, 7.197, 4.983, 6.423, 4.394, 5.914, 4.03, 5.598, 3.79, 5.411, 3.595, 5.225) T[5, .] = (50, 6.74, 6.74, 4.19, 4.94, 3.333, 4.313, 2.873, 3.973, 2.614, 3.746, 2.435, 3.6, 2.309, 3.507, 2.205, 3.421, 8.51, 8.51, 5.22, 6.07, 4.07, 5.19, 3.5, 4.7, 3.136, 4.416, 2.9, 4.218, 2.726, 4.057, 2.593, 3.941, 12.73, 12.73, 7.56, 8.685, 5.817, 7.303, 4.865, 6.36, 4.306, 5.874, 3.955, 5.583, 3.656, 5.331, 3.498, 5.149) T[6, .] = (55, 6.7, 6.7, 4.155, 4.925, 3.28, 4.273, 2.843, 3.92, 2.578, 3.71, 2.393, 3.583, 2.27, 3.486, 2.181, 3.398, 8.39, 8.39, 5.125, 6.045, 3.987, 5.09, 3.408, 4.623, 3.068, 4.334, 2.848, 4.16, 2.676, 3.999, 2.556, 3.904, 12.7, 12.7, 7.435, 8.46, 5.707, 6.977, 4.828, 6.195, 4.244, 5.726, 3.928, 5.408, 3.636, 5.169, 3.424, 4.989) T[7, .] = (60, 6.7, 6.7, 4.145, 4.95, 3.27, 4.26, 2.838, 3.923, 2.568, 3.712, 2.385, 3.565, 2.253, 3.436, 2.155, 3.353, 8.46, 8.46, 5.125, 6, 4, 5.057, 3.415, 4.615, 3.062, 4.314, 2.817, 4.097, 2.643, 3.939, 2.513, 3.823, 12.49, 12.49, 7.4, 8.51, 5.697, 6.987, 4.748, 6.188, 4.176, 5.676, 3.783, 5.338, 3.531, 5.081, 3.346, 4.895) T[8, .] = (65, 6.74, 6.74, 4.175, 4.93, 3.3, 4.25, 2.843, 3.923, 2.574, 3.682, 2.397, 3.543, 2.256, 3.43, 2.156, 3.334, 8.49, 8.49, 5.13, 5.98, 4.01, 5.08, 3.435, 4.583, 3.068, 4.274, 2.835, 4.09, 2.647, 3.921, 2.525, 3.808, 12.4, 12.4, 7.32, 8.435, 5.583, 6.853, 4.69, 6.143, 4.188, 5.694, 3.783, 5.3, 3.501, 5.051, 3.31, 4.871) T[9, .] = (70, 6.67, 6.67, 4.125, 4.88, 3.25, 4.237, 2.818, 3.88, 2.552, 3.648, 2.363, 3.51, 2.233, 3.407, 2.138, 3.325, 8.37, 8.37, 5.055, 5.915, 3.947, 5.02, 3.37, 4.545, 3.022, 4.256, 2.788, 4.073, 2.629, 3.906, 2.494, 3.786, 12.24, 12.24, 7.17, 8.405, 5.487, 6.88, 4.635, 6.055, 4.098, 5.57, 3.747, 5.285, 3.436, 5.044, 3.261, 4.821) T[10, .] = (75, 6.72, 6.72, 4.15, 4.885, 3.277, 4.243, 2.838, 3.898, 2.558, 3.654, 2.38, 3.515, 2.244, 3.397, 2.134, 3.313, 8.42, 8.42, 5.14, 5.92, 3.983, 5.06, 3.408, 4.55, 3.042, 4.244, 2.802, 4.065, 2.637, 3.9, 2.503, 3.768, 12.54, 12.54, 7.225, 8.3, 5.513, 6.86, 4.725, 6.08, 4.168, 5.548, 3.772, 5.213, 3.496, 4.966, 3.266, 4.801) T[11, .] = (80, 6.72, 6.72, 4.135, 4.895, 3.26, 4.247, 2.823, 3.885, 2.548, 3.644, 2.355, 3.5, 2.236, 3.381, 2.129, 3.289, 8.4, 8.4, 5.06, 5.93, 3.94, 5.043, 3.363, 4.515, 3.01, 4.216, 2.787, 4.015, 2.627, 3.864, 2.476, 3.746, 12.12, 12.12, 7.095, 8.26, 5.407, 6.783, 4.568, 5.96, 4.096, 5.512, 3.725, 5.163, 3.457, 4.943, 3.233, 4.76) } if (cs == 4) { T = J(10, 49, .) T[1, .] = (30, 5.785, 5.785, 4.427, 4.957, 3.77, 4.535, 3.378, 4.274, 3.097, 4.118, 2.907, 4.01, 2.781, 3.941, 2.681, 3.887, 7.04, 7.04, 5.377, 5.963, 4.535, 5.415, 4.048, 5.09, 3.715, 4.878, 3.504, 4.743, 3.326, 4.653, 3.194, 4.604, 10.2, 10.2, 7.593, 8.35, 6.428, 7.505, 5.666, 6.988, 5.205, 6.64, 4.85, 6.473, 4.689, 6.358, 4.49, 6.328) T[2, .] = (35, 5.69, 5.69, 4.38, 4.867, 3.698, 4.42, 3.29, 4.176, 3.035, 3.997, 2.831, 3.879, 2.685, 3.785, 2.578, 3.71, 6.9, 6.9, 5.233, 5.777, 4.433, 5.245, 3.936, 4.918, 3.578, 4.668, 3.353, 4.5, 3.174, 4.383, 3.057, 4.319, 9.975, 9.975, 7.477, 8.213, 6.328, 7.408, 5.654, 6.926, 5.147, 6.617, 4.849, 6.511, 4.629, 5.698, 4.489, 5.064) T[3, .] = (40, 5.68, 5.68, 4.343, 4.823, 3.663, 4.378, 3.264, 4.094, 2.985, 3.918, 2.781, 3.813, 2.634, 3.719, 2.517, 3.65, 6.87, 6.87, 5.18, 5.733, 4.36, 5.138, 3.85, 4.782, 3.512, 4.587, 3.257, 4.431, 3.07, 4.309, 2.933, 4.224, 9.575, 9.575, 7.207, 7.86, 5.98, 6.973, 5.258, 6.526, 4.763, 6.2, 4.427, 5.837, 4.154, 5.699, 3.971, 5.486) T[4, .] = (45, 5.625, 5.625, 4.3, 4.78, 3.625, 4.33, 3.226, 4.054, 2.95, 3.862, 2.75, 3.739, 2.606, 3.644, 2.484, 3.57, 6.75, 6.75, 5.13, 5.68, 4.335, 5.078, 3.822, 4.714, 3.47, 4.47, 3.211, 4.309, 3.025, 4.198, 2.899, 4.087, 9.555, 9.555, 7.133, 7.82, 5.878, 6.87, 5.15, 6.28, 4.628, 5.865, 4.251, 5.596, 3.998, 5.463, 3.829, 5.313) T[5, .] = (50, 5.57, 5.57, 4.23, 4.74, 3.573, 4.288, 3.174, 4.004, 2.905, 3.822, 2.703, 3.697, 2.55, 3.609, 2.44, 3.523, 6.685, 6.685, 5.043, 5.607, 4.225, 5.03, 3.73, 4.666, 3.383, 4.432, 3.149, 4.293, 2.975, 4.143, 2.832, 4.012, 9.32, 9.32, 7.017, 7.727, 5.805, 6.79, 5.05, 6.182, 4.557, 5.793, 4.214, 5.52, 3.983, 5.345, 3.762, 5.172) T[6, .] = (55, 5.57, 5.57, 4.23, 4.73, 3.553, 4.238, 3.132, 3.956, 2.868, 3.782, 2.674, 3.659, 2.538, 3.56, 2.42, 3.481, 6.66, 6.66, 5.013, 5.547, 4.183, 4.955, 3.692, 4.582, 3.358, 4.365, 3.131, 4.206, 2.946, 4.065, 2.791, 3.95, 9.3, 9.3, 6.893, 7.537, 5.678, 6.578, 4.99, 6.018, 4.455, 5.615, 4.111, 5.329, 3.87, 5.171, 3.643, 5.021) T[7, .] = (60, 5.555, 5.555, 4.203, 4.693, 3.54, 4.235, 3.13, 3.968, 2.852, 3.773, 2.653, 3.637, 2.51, 3.519, 2.392, 3.444, 6.63, 6.63, 4.98, 5.527, 4.18, 4.938, 3.684, 4.584, 3.323, 4.333, 3.086, 4.154, 2.9, 3.999, 2.756, 3.892, 9.245, 9.245, 6.78, 7.377, 5.62, 6.503, 4.928, 5.95, 4.412, 5.545, 4.013, 5.269, 3.775, 5.086, 3.584, 4.922) T[8, .] = (65, 5.51, 5.51, 4.187, 4.66, 3.535, 4.208, 3.122, 3.942, 2.848, 3.743, 2.647, 3.603, 2.499, 3.49, 2.379, 3.406, 6.55, 6.55, 4.95, 5.467, 4.123, 4.903, 3.626, 4.538, 3.3, 4.28, 3.063, 4.123, 2.88, 3.978, 2.73, 3.879, 8.96, 8.96, 6.707, 7.36, 5.545, 6.453, 4.848, 5.842, 4.347, 5.552, 4.02, 5.263, 3.758, 5.04, 3.557, 4.902) T[9, .] = (70, 5.53, 5.53, 4.173, 4.647, 3.505, 4.198, 3.098, 3.92, 2.832, 3.738, 2.631, 3.589, 2.485, 3.473, 2.363, 3.394, 6.53, 6.53, 4.93, 5.457, 4.1, 4.9, 3.6, 4.512, 3.272, 4.272, 3.043, 4.1, 2.86, 3.951, 2.711, 3.842, 8.89, 8.89, 6.577, 7.313, 5.448, 6.435, 4.76, 5.798, 4.293, 5.46, 3.966, 5.234, 3.72, 5.004, 3.509, 4.808) T[10, .] = (75, 5.53, 5.53, 4.193, 4.647, 3.505, 4.213, 3.11, 3.9, 2.832, 3.717, 2.636, 3.579, 2.486, 3.469, 2.372, 3.37, 6.58, 6.58, 4.937, 5.443, 4.12, 4.855, 3.624, 4.488, 3.298, 4.26, 3.054, 4.079, 2.874, 3.914, 2.718, 3.807, 8.905, 8.905, 6.613, 7.253, 5.505, 6.298, 4.808, 5.786, 4.3, 5.377, 3.984, 5.153, 3.728, 4.954, 3.511, 4.789) } if (cs == 5) { T = J(11, 49, .) T[1, .] = (30, 10.34, 10.34, 6.01, 6.78, 4.577, 5.6, 3.868, 4.965, 3.43, 4.624, 3.157, 4.412, 2.977, 4.26, 2.843, 4.16, 12.74, 12.74, 7.36, 8.265, 5.55, 6.747, 4.683, 5.98, 4.154, 5.54, 3.818, 5.253, 3.576, 5.065, 3.394, 4.939, 18.56, 18.56, 10.605, 11.65, 7.977, 9.413, 6.643, 8.313, 5.856, 7.578, 5.347, 7.242, 5.046, 6.93, 4.779, 6.821) T[2, .] = (35, 10.24, 10.24, 5.95, 6.68, 4.517, 5.48, 3.8, 4.888, 3.374, 4.512, 3.087, 4.277, 2.879, 4.114, 2.729, 3.985, 12.58, 12.58, 7.21, 8.055, 5.457, 6.57, 4.568, 5.795, 4.036, 5.304, 3.673, 5.002, 3.426, 4.79, 3.251, 4.64, 18.02, 18.02, 10.365, 11.295, 7.643, 9.063, 6.38, 7.73, 5.604, 7.172, 5.095, 6.77, 4.704, 6.537, 4.459, 6.206) T[3, .] = (40, 10.16, 10.16, 5.915, 6.63, 4.477, 5.42, 3.76, 4.795, 3.334, 4.438, 3.032, 4.213, 2.831, 4.04, 2.668, 3.92, 12.51, 12.51, 7.135, 7.98, 5.387, 6.437, 4.51, 5.643, 3.958, 5.226, 3.577, 4.923, 3.327, 4.7, 3.121, 4.564, 17.91, 17.91, 10.15, 11.23, 7.527, 8.803, 6.238, 7.74, 5.376, 7.092, 4.885, 6.55, 4.527, 6.263, 4.31, 5.965) T[4, .] = (45, 10.15, 10.15, 5.88, 6.64, 4.437, 5.377, 3.74, 4.78, 3.298, 4.378, 3.012, 4.147, 2.796, 3.97, 2.635, 3.838, 12.4, 12.4, 7.08, 7.91, 5.36, 6.373, 4.45, 5.56, 3.89, 5.104, 3.532, 4.8, 3.267, 4.584, 3.091, 4.413, 17.5, 17.5, 9.89, 10.965, 7.317, 8.72, 6.053, 7.458, 5.224, 6.696, 4.715, 6.293, 4.364, 6.006, 4.109, 5.785) T[5, .] = (50, 10.02, 10.02, 5.78, 6.54, 4.38, 5.35, 3.673, 4.715, 3.24, 4.35, 2.95, 4.11, 2.75, 3.944, 2.59, 3.789, 12.17, 12.17, 6.985, 7.86, 5.247, 6.303, 4.368, 5.545, 3.834, 5.064, 3.48, 4.782, 3.229, 4.536, 3.039, 4.339, 17.53, 17.53, 9.895, 10.965, 7.337, 8.643, 5.995, 7.335, 5.184, 6.684, 4.672, 6.232, 4.31, 5.881, 4.055, 5.64) T[6, .] = (55, 10.11, 10.11, 5.8, 6.515, 4.37, 5.303, 3.64, 4.67, 3.21, 4.294, 2.927, 4.068, 2.724, 3.893, 2.573, 3.76, 12.17, 12.17, 6.93, 7.785, 5.19, 6.223, 4.313, 5.425, 3.794, 4.986, 3.442, 4.69, 3.197, 4.46, 2.989, 4.271, 17.48, 17.48, 9.8, 10.675, 7.227, 8.34, 5.955, 7.225, 5.108, 6.494, 4.608, 5.977, 4.23, 5.713, 3.955, 5.474) T[7, .] = (60, 10.03, 10.03, 5.765, 6.5, 4.35, 5.283, 3.645, 4.678, 3.2, 4.31, 2.912, 4.047, 2.709, 3.856, 2.551, 3.716, 12.2, 12.2, 6.905, 7.735, 5.19, 6.2, 4.298, 5.445, 3.772, 4.956, 3.407, 4.632, 3.137, 4.393, 2.956, 4.23, 17.02, 17.02, 9.585, 10.42, 7.057, 8.243, 5.835, 7.108, 5.066, 6.394, 4.505, 5.92, 4.117, 5.597, 3.87, 5.338) T[8, .] = (65, 9.97, 9.97, 5.755, 6.47, 4.353, 5.257, 3.638, 4.643, 3.196, 4.262, 2.897, 4.022, 2.69, 3.83, 2.531, 3.685, 11.96, 11.96, 6.89, 7.66, 5.137, 6.173, 4.268, 5.415, 3.732, 4.92, 3.372, 4.613, 3.137, 4.363, 2.924, 4.206, 16.85, 16.85, 9.475, 10.515, 7.013, 8.23, 5.795, 7.053, 4.974, 6.378, 4.482, 5.923, 4.111, 5.586, 3.835, 5.339) T[9, .] = (70, 10.02, 10.02, 5.765, 6.455, 4.33, 5.243, 3.615, 4.635, 3.182, 4.258, 2.893, 4.008, 2.683, 3.807, 2.519, 3.669, 12, 12, 6.86, 7.645, 5.11, 6.19, 4.235, 5.363, 3.72, 4.904, 3.368, 4.59, 3.107, 4.343, 2.913, 4.168, 16.66, 16.66, 9.37, 10.32, 6.873, 8.163, 5.663, 6.953, 4.922, 6.328, 4.428, 5.898, 4.07, 5.534, 3.774, 5.248) T[10, .] = (75, 10.03, 10.03, 5.765, 6.47, 4.323, 5.273, 3.618, 4.63, 3.182, 4.248, 2.89, 3.993, 2.681, 3.8, 2.53, 3.648, 12.08, 12.08, 6.88, 7.675, 5.14, 6.153, 4.253, 5.333, 3.724, 4.88, 3.382, 4.567, 3.111, 4.31, 2.915, 4.143, 16.61, 16.61, 9.325, 10.325, 6.93, 8.027, 5.698, 6.97, 4.932, 6.224, 4.393, 5.788, 4.06, 5.459, 3.768, 5.229) T[11, .] = (80, 9.96, 9.96, 5.725, 6.45, 4.307, 5.223, 3.588, 4.605, 3.16, 4.23, 2.867, 3.975, 2.657, 3.776, 2.504, 3.631, 12.06, 12.06, 6.82, 7.67, 5.067, 6.103, 4.203, 5.32, 3.678, 4.84, 3.335, 4.535, 3.077, 4.284, 2.885, 4.111, 16.6, 16.6, 9.17, 10.24, 6.73, 8.053, 5.62, 6.908, 4.89, 6.164, 4.375, 5.703, 4, 5.397, 3.728, 5.16) } nr = rows(T) lo = 1 hi = 1 w = 0 if (n > T[1, 1]) { if (n >= T[nr, 1]) { lo = nr hi = nr w = 0 } else { for (i = 1; i <= nr - 1; i++) { if (n >= T[i, 1]) { lo = i } } hi = lo + 1 w = (n - T[lo, 1]) / (T[hi, 1] - T[lo, 1]) } } a = T[lo, .] b = T[hi, .] O = J(3, 2, .) for (i = 1; i <= 3; i++) { j = 1 + (i - 1) * 16 + k * 2 O[i, 1] = (1 - w) * a[j + 1] + w * b[j + 1] O[i, 2] = (1 - w) * a[j + 2] + w * b[j + 2] } return(O) } // ------------------------------------------------------------------ // Bounds decision code, equation (118): // 2 = reject H0 (statistic beyond the I(1) bound) // 1 = inconclusive // 0 = do not reject (statistic inside the I(0) bound) // . = bounds unavailable // upper != 0 for the right-tailed F statistic; upper == 0 for the // left-tailed t statistic on lambda. // ------------------------------------------------------------------ real scalar _surardl_decide(real scalar stat, real scalar c0, real scalar c1, real scalar upper) { if (c0 >= . | c1 >= .) { return(.) } if (upper != 0) { if (stat > c1) { return(2) } if (stat < c0) { return(0) } return(1) } if (stat < c1) { return(2) } if (stat > c0) { return(0) } return(1) } // ================================================================ // from _sa_engine.mata // ================================================================ real scalar _surardl_ver() { return(101) } struct surardl_eq { real colvector y real matrix X real scalar m real scalar p real rowvector q real rowvector rst real scalar ilam real rowvector ibeta real scalar itr real scalar ico real rowvector ipsi real rowvector igam real scalar off string rowvector cn } // Column of unit i, variable v (v = 0 is the dependent variable, // v = 1..K the regressors) in a level-series matrix LV. real scalar _surardl_col(real scalar i, real scalar v, real scalar K) { return((i - 1) * (K + 1) + 1 + v) } // Number of initial observations consumed by the lag structure. real scalar _surardl_lagmax(real matrix LAGS, real scalar K) { real scalar i, k, L L = 1 for (i = 1; i <= rows(LAGS); i++) { if (LAGS[i, 1] > L) { L = LAGS[i, 1] } for (k = 1; k <= K; k++) { if (LAGS[i, k+1] > L) { L = LAGS[i, k+1] } } } return(L) } // ------------------------------------------------------------------ // Build the per-equation conditional ECM designs directly from the // matrix of level series LV (Tfull x N*(K+1), gap-free and balanced; // column block per unit, ordered y_i, x_1i, ..., x_Ki) and the vector // of trend values TR (Tfull x 1). // Implements equations (7) and (12)-(15) of Arif & Furuoka (2026). // The same routine serves the observed data and every bootstrap // replication, so the bootstrap reproduces the whole statistic. // ------------------------------------------------------------------ struct surardl_eq vector _surardl_build(real matrix LV, real colvector TR, real scalar N, real scalar K, real matrix LAGS, real scalar cs, string scalar dep, string rowvector xn) { struct surardl_eq vector E real scalar i, k, j, off, pos, L, T, Tf, c0, ck real colvector rr string rowvector cn E = surardl_eq(N) Tf = rows(LV) L = _surardl_lagmax(LAGS, K) T = Tf - L rr = ((L + 1) :: Tf) off = 0 for (i = 1; i <= N; i++) { E[i].p = LAGS[i, 1] E[i].q = LAGS[i, (2..(K+1))] c0 = _surardl_col(i, 0, K) E[i].y = LV[rr, c0] - LV[rr :- 1, c0] // --- lagged levels: lambda then beta_1..beta_K (eq. 14) --- E[i].X = LV[rr :- 1, c0] cn = J(1, 1, "L." + dep) for (k = 1; k <= K; k++) { ck = _surardl_col(i, k, K) E[i].X = E[i].X , LV[rr :- 1, ck] cn = cn , ("L." + xn[k]) } E[i].ilam = 1 E[i].ibeta = (2..(K+1)) pos = K + 1 // --- lagged differences of the dependent variable (psi, eq. 10) --- E[i].ipsi = J(1, 0, .) for (j = 1; j <= E[i].p - 1; j++) { E[i].X = E[i].X , (LV[rr :- j, c0] - LV[rr :- (j+1), c0]) pos = pos + 1 E[i].ipsi = E[i].ipsi , pos if (j == 1) { cn = cn , ("LD." + dep) } else { cn = cn , ("L" + strofreal(j) + "D." + dep) } } // --- current and lagged differences of regressors (gamma, eq. 11) --- E[i].igam = J(1, 0, .) for (k = 1; k <= K; k++) { ck = _surardl_col(i, k, K) for (j = 0; j <= E[i].q[k] - 1; j++) { E[i].X = E[i].X , (LV[rr :- j, ck] - LV[rr :- (j+1), ck]) pos = pos + 1 E[i].igam = E[i].igam , pos if (j == 0) { cn = cn , ("D." + xn[k]) } else { if (j == 1) { cn = cn , ("LD." + xn[k]) } else { cn = cn , ("L" + strofreal(j) + "D." + xn[k]) } } } } // --- deterministic terms --- E[i].itr = 0 E[i].ico = 0 if (cs >= 4) { E[i].X = E[i].X , TR[rr] pos = pos + 1 E[i].itr = pos cn = cn , "_trend" } if (cs >= 2) { E[i].X = E[i].X , J(T, 1, 1) pos = pos + 1 E[i].ico = pos cn = cn , "_cons" } E[i].m = pos E[i].cn = cn E[i].off = off off = off + pos // --- restriction set for H0: lambda_i = beta_1i = ... = beta_Ki = 0, // plus the restricted deterministic term in cases 2 and 4. // eq. (61)/(63); see also Kripfganz & Schneider (2020, fn. 12): // q = K+1 in cases 1, 3, 5 and q = K+2 in cases 2, 4. E[i].rst = (1..(K+1)) if (cs == 2) { E[i].rst = E[i].rst , E[i].ico } if (cs == 4) { E[i].rst = E[i].rst , E[i].itr } } return(E) } real scalar _surardl_M(struct surardl_eq vector E) { real scalar i, M M = 0 for (i = 1; i <= length(E); i++) { M = M + E[i].m } return(M) } // ------------------------------------------------------------------ // Equation-by-equation OLS, equation (18). Returns the T x N residual // matrix and fills bols, seols, rssv, r2v, rmsev. // ------------------------------------------------------------------ real matrix _surardl_ols(struct surardl_eq vector E, real colvector bols, real colvector seols, real colvector rssv, real colvector r2v, real colvector rmsev) { real scalar i, N, T, s2, tss real matrix U, iXX real colvector bi, ui, my N = length(E) T = rows(E[1].y) U = J(T, N, 0) bols = J(0, 1, .) seols = J(0, 1, .) rssv = J(N, 1, .) r2v = J(N, 1, .) rmsev = J(N, 1, .) for (i = 1; i <= N; i++) { iXX = invsym(quadcross(E[i].X, E[i].X)) bi = iXX * quadcross(E[i].X, E[i].y) ui = E[i].y - E[i].X * bi U[., i] = ui rssv[i] = quadcross(ui, ui) s2 = rssv[i] / (T - E[i].m) rmsev[i]= sqrt(s2) bols = bols \ bi seols = seols \ sqrt(diagonal(s2 * iXX)) my = E[i].y :- mean(E[i].y) tss = quadcross(my, my) r2v[i] = 1 - rssv[i] / tss } return(U) } // ------------------------------------------------------------------ // Contemporaneous covariance matrix Sigma-hat, equations (55)-(56). // dfk = 0 : divide by T (Zellner MLE; sureg default; the // divisor used by the authors' code) // dfk = 1 : T - max(m_i, m_j) (paper's third option) // dfk = 2 : sqrt((T-m_i)(T-m_j)) (sureg's dfk) // dfk = 3 : T - m_i - m_j (sureg's dfk2, leading term) // ------------------------------------------------------------------ real matrix _surardl_sigma(real matrix U, struct surardl_eq vector E, real scalar dfk) { real scalar i, j, N, T, d real matrix S N = cols(U) T = rows(U) S = J(N, N, 0) for (i = 1; i <= N; i++) { for (j = 1; j <= N; j++) { d = T if (dfk == 1) { d = T - max((E[i].m, E[j].m)) } if (dfk == 2) { d = sqrt((T - E[i].m) * (T - E[j].m)) } if (dfk == 3) { d = T - E[i].m - E[j].m if (d <= 0) { d = T } } S[i, j] = quadcross(U[., i], U[., j]) / d } } return(S) } // ------------------------------------------------------------------ // FGLS cross-products, equations (50)-(53). The (i,j) block of A is // s^ij X_i'X_j and the i-th block of b is sum_j s^ij X_i'y_j. // ------------------------------------------------------------------ void _surardl_Ab(struct surardl_eq vector E, real matrix Si, real matrix A, real colvector b) { real scalar i, j, N, M, r0, c0 real colvector bi N = length(E) M = _surardl_M(E) A = J(M, M, 0) b = J(M, 1, 0) for (i = 1; i <= N; i++) { r0 = E[i].off bi = J(E[i].m, 1, 0) for (j = 1; j <= N; j++) { c0 = E[j].off A[|r0+1, c0+1 \ r0+E[i].m, c0+E[j].m|] = Si[i, j] * quadcross(E[i].X, E[j].X) bi = bi + Si[i, j] * quadcross(E[i].X, E[j].y) } b[|r0+1 \ r0+E[i].m|] = bi } } void _surardl_resid(struct surardl_eq vector E, real colvector th, real matrix U) { real scalar i, N N = length(E) for (i = 1; i <= N; i++) { U[., i] = E[i].y - E[i].X * th[|E[i].off+1 \ E[i].off+E[i].m|] } } // ------------------------------------------------------------------ // Two-step / iterated FGLS, equations (57)-(58) and the iteration // scheme (a)-(d) of section 2.3. // ------------------------------------------------------------------ void _surardl_fgls(struct surardl_eq vector E, real matrix U0, real scalar dfk, real scalar iterate, real scalar tol, real colvector th, real matrix V, real matrix S, real matrix U, real scalar nit, real scalar conv) { real scalar it, dmax real matrix A, Si real colvector b, thold U = U0 S = _surardl_sigma(U, E, dfk) Si = invsym(S) _surardl_Ab(E, Si, A, b) V = invsym(A) th = V * b _surardl_resid(E, th, U) nit = 1 conv = 1 if (iterate > 0) { conv = 0 for (it = 1; it <= iterate; it++) { thold = th S = _surardl_sigma(U, E, dfk) Si = invsym(S) _surardl_Ab(E, Si, A, b) V = invsym(A) th = V * b _surardl_resid(E, th, U) nit = it dmax = max(abs(th - thold)) if (dmax < tol) { conv = 1 it = iterate + 1 } } } } // ------------------------------------------------------------------ // System-level restriction matrix for H0^(i), equations (63)-(64). // ------------------------------------------------------------------ real matrix _surardl_R(struct surardl_eq vector E, real scalar i, real scalar M) { real scalar r, nr real matrix R nr = length(E[i].rst) R = J(nr, M, 0) for (r = 1; r <= nr; r++) { R[r, E[i].off + E[i].rst[r]] = 1 } return(R) } // Restriction matrix selecting only the lagged x levels (beta_1..beta_K) // of equation i. This is the second PSS degenerate case, tested by // McNown, Sam & Goh (2018) with t_IDV. real matrix _surardl_Rbeta(struct surardl_eq vector E, real scalar i, real scalar M) { real scalar r, nr real matrix R nr = length(E[i].ibeta) R = J(nr, M, 0) for (r = 1; r <= nr; r++) { R[r, E[i].off + E[i].ibeta[r]] = 1 } return(R) } // ------------------------------------------------------------------ // Wald statistic, equations (77) / (112); F form F = W/q, eq. (76). // Returns (W, q, F). // ------------------------------------------------------------------ real rowvector _surardl_wald(real matrix R, real colvector th, real matrix V) { real colvector d real matrix Q real scalar w, q d = R * th Q = R * V * R' w = (d' * invsym(Q) * d) q = rows(R) return((w, q, w / q)) } // ------------------------------------------------------------------ // RSS-based F statistic, equations (73)-(76). The restricted system // drops the level terms from equation i only (eq. 69: "the models for // all other units j != i remain unrestricted"). With allq != 0 the // level terms are dropped from every equation simultaneously, which is // what the authors' own OxGauss / Stata scripts do. // Sigma-hat is held at its unrestricted value, as in eq. (73). // ------------------------------------------------------------------ real colvector _surardl_rssF(struct surardl_eq vector E, real matrix Si, real colvector rssU, real scalar allq) { real scalar i, j, N, T, M, nk, q real colvector F, bR, keep, rssR, b, uj real matrix A, Si2 struct surardl_eq vector ER N = length(E) T = rows(E[1].y) F = J(N, 1, .) for (i = 1; i <= N; i++) { ER = E for (j = 1; j <= N; j++) { if (allq != 0 | j == i) { keep = J(0, 1, .) for (nk = 1; nk <= E[j].m; nk++) { if (!anyof(E[j].rst, nk)) { keep = keep \ nk } } ER[j].X = E[j].X[., keep] ER[j].m = rows(keep) } } // recompute offsets M = 0 for (j = 1; j <= N; j++) { ER[j].off = M M = M + ER[j].m } Si2 = Si _surardl_Ab(ER, Si2, A, b) bR = invsym(A) * b rssR = J(N, 1, .) for (j = 1; j <= N; j++) { uj = ER[j].y - ER[j].X * bR[|ER[j].off+1 \ ER[j].off+ER[j].m|] rssR[j] = quadcross(uj, uj) } q = length(E[i].rst) F[i] = ((rssR[i] - rssU[i]) / q) / (rssU[i] / (T - E[i].m)) } return(F) } // ------------------------------------------------------------------ // Long-run coefficients, equation (84), with delta-method variances, // equations (85)-(91). // Returns a (N*K) x 5 matrix: unit, regressor, theta, se, and the // covariance block determinant flag. // ------------------------------------------------------------------ real matrix _surardl_longrun(struct surardl_eq vector E, real colvector th, real matrix V, real scalar K) { real scalar i, k, N, il, ib, lam, bet, vll, vbb, vlb, va real matrix Out N = length(E) Out = J(N * K, 4, .) for (i = 1; i <= N; i++) { il = E[i].off + E[i].ilam lam = th[il] vll = V[il, il] for (k = 1; k <= K; k++) { ib = E[i].off + E[i].ibeta[k] bet = th[ib] vbb = V[ib, ib] vlb = V[il, ib] // eq. (89): Var = (b^2/l^4) v_ll - (2b/l^3) v_lb + (1/l^2) v_bb va = (bet * bet) / (lam^4) * vll va = va - 2 * bet / (lam^3) * vlb va = va + vbb / (lam * lam) Out[(i-1)*K + k, 1] = i Out[(i-1)*K + k, 2] = k Out[(i-1)*K + k, 3] = -bet / lam Out[(i-1)*K + k, 4] = sqrt(va) } } return(Out) } // ------------------------------------------------------------------ // Fieller (1954) confidence set for theta = -beta/lambda. Inverts // H0: beta + theta*lambda = 0 and therefore never divides by lambda. // Recommended over the delta method in section 4.6 of the paper. // Returns (lo, hi, type) with type // 0 = bounded interval // 1 = unbounded: complement of (lo,hi), i.e. (-inf,lo] U [hi,+inf) // 2 = entire real line // 3 = empty set // ------------------------------------------------------------------ real rowvector _surardl_fieller(real scalar lam, real scalar bet, real scalar vll, real scalar vbb, real scalar vlb, real scalar z) { real scalar A, B, C, Disc, r1, r2, zz zz = z * z A = lam * lam - zz * vll B = 2 * (bet * lam - zz * vlb) C = bet * bet - zz * vbb Disc = B * B - 4 * A * C if (A > 0) { if (Disc <= 0) { return((., ., 3)) } r1 = (-B - sqrt(Disc)) / (2 * A) r2 = (-B + sqrt(Disc)) / (2 * A) return((min((r1, r2)), max((r1, r2)), 0)) } if (A < 0) { if (Disc <= 0) { return((., ., 2)) } r1 = (-B - sqrt(Disc)) / (2 * A) r2 = (-B + sqrt(Disc)) / (2 * A) return((min((r1, r2)), max((r1, r2)), 1)) } return((., ., 2)) } // ------------------------------------------------------------------ // Breusch & Pagan (1980) LM test of diagonal Sigma, equations (96)-(99), // and the Pesaran (2004) CD statistic, equation (100). // NOTE on eq. (100): as printed in the draft the CD statistic omits the // factor 2 inside the square root. The formula implemented here is // Pesaran's own, CD = sqrt(2T/(N(N-1))) * sum_{i= T - 2) { return(Out) } for (i = 1; i <= N; i++) { u = U[((L+1)::T), i] W = E[i].X[((L+1)::T), .] for (l = 1; l <= L; l++) { ui = U[((L+1-l)::(T-l)), i] W = W , ui } n = rows(u) e2 = u - W * (invsym(quadcross(W, W)) * quadcross(W, u)) mu = u :- mean(u) tss = quadcross(mu, mu) r2 = 1 - quadcross(e2, e2) / tss Out[i, 1] = n * r2 Out[i, 2] = L Out[i, 3] = chi2tail(L, n * r2) } return(Out) } // ================================================================ // from _sa_boot.mata // ================================================================ // ================================================================== // Statistic vector for one data set. Returns N x 3: // col 1 F_i bounds F for H0: lambda_i = beta_.i = 0 (plus the // restricted deterministic term in cases 2 and 4), // equations (76)-(77) // col 2 t_i t statistic on lambda_i (PSS step 2; first // degenerate case) // col 3 Fidv_i F for H0: beta_1i = ... = beta_Ki = 0 (second // degenerate case; t_IDV of McNown et al. 2018) // ================================================================== real matrix _surardl_stats(struct surardl_eq vector E, real scalar dfk, real scalar iterate, real scalar tol) { real scalar i, N, M, il, nit, conv real matrix S, U, U0, V, Out, R real colvector th, bols, seols, rssv, r2v, rmsev real rowvector w N = length(E) M = _surardl_M(E) U0 = _surardl_ols(E, bols, seols, rssv, r2v, rmsev) _surardl_fgls(E, U0, dfk, iterate, tol, th, V, S, U, nit, conv) Out = J(N, 3, .) for (i = 1; i <= N; i++) { R = _surardl_R(E, i, M) w = _surardl_wald(R, th, V) Out[i, 1] = w[3] il = E[i].off + E[i].ilam Out[i, 2] = th[il] / sqrt(V[il, il]) R = _surardl_Rbeta(E, i, M) w = _surardl_wald(R, th, V) Out[i, 3] = w[3] } return(Out) } // Trap a failed replication (singular design or non-invertible // Sigma-hat) instead of aborting the whole bootstrap. real scalar _surardl_try(struct surardl_eq vector E, real scalar dfk, real scalar iterate, real scalar tol, real matrix St) { real scalar i, N, ok N = length(E) ok = 1 for (i = 1; i <= N; i++) { if (rank(quadcross(E[i].X, E[i].X)) < E[i].m) { ok = 0 } } if (ok == 0) { return(1) } St = _surardl_stats(E, dfk, iterate, tol) if (hasmissing(St)) { return(1) } return(0) } // Rescale and recentre residuals column by column; McNown, Sam & Goh // (2018, eq. 13), following Davidson & MacKinnon (2004). np holds the // number of estimated parameters behind each column. real matrix _surardl_rescale(real matrix U, real rowvector np) { real scalar n, c, f real matrix R n = rows(U) R = U for (c = 1; c <= cols(U); c++) { f = 1 if (n > np[c]) { f = sqrt(n / (n - np[c])) } R[., c] = f :* (U[., c] :- mean(U[., c])) } return(R) } // ================================================================== // Marginal (conditioning) model for the regressors of unit i, used // only to generate pseudo-data. Equation X of McNown, Sam & Goh // (2018, eq. 12), generalised to K regressors and VECM order mo: // // dx_{k,i,t} = a + b'(y_{i,t-1}, x_{i,t-1}') // + sum_{j=1..mo} phi_j'(dy_{i,t-j}, dx_{i,t-j}') // + v_{k,i,t} // // Fitted by OLS on rows Lg+1..Tf. Returns the (nreg x K) coefficient // matrix and fills VR with the (Tf-Lg) x K residuals. // ================================================================== real matrix _surardl_marg(real matrix LV, real colvector TR, real scalar i, real scalar K, real scalar Lg, real scalar mo, real scalar cs, real matrix VR) { real scalar k, j, Tb, Tf, c0, ck real colvector rr, yk real matrix W, B Tf = rows(LV) Tb = Tf - Lg rr = ((Lg + 1) :: Tf) c0 = _surardl_col(i, 0, K) W = LV[rr :- 1, c0] for (k = 1; k <= K; k++) { ck = _surardl_col(i, k, K) W = W , LV[rr :- 1, ck] } for (j = 1; j <= mo; j++) { W = W , (LV[rr :- j, c0] - LV[rr :- (j+1), c0]) for (k = 1; k <= K; k++) { ck = _surardl_col(i, k, K) W = W , (LV[rr :- j, ck] - LV[rr :- (j+1), ck]) } } if (cs >= 4) { W = W , TR[rr] } if (cs >= 2) { W = W , J(Tb, 1, 1) } B = J(cols(W), K, 0) VR = J(Tb, K, 0) for (k = 1; k <= K; k++) { ck = _surardl_col(i, k, K) yk = LV[rr, ck] - LV[rr :- 1, ck] B[., k] = invsym(quadcross(W, W)) * quadcross(W, yk) VR[., k] = yk - W * B[., k] } return(B) } // Regressor row for the marginal model of unit i at time t, built // from the pseudo-level matrix LVs. Must match _surardl_marg()'s W. real colvector _surardl_margrow(real matrix LVs, real colvector TR, real scalar i, real scalar K, real scalar t, real scalar mo, real scalar cs) { real scalar k, j, c0 real colvector z c0 = _surardl_col(i, 0, K) z = LVs[t-1, c0] for (k = 1; k <= K; k++) { z = z \ LVs[t-1, _surardl_col(i, k, K)] } for (j = 1; j <= mo; j++) { z = z \ (LVs[t-j, c0] - LVs[t-j-1, c0]) for (k = 1; k <= K; k++) { z = z \ (LVs[t-j, _surardl_col(i, k, K)] - LVs[t-j-1, _surardl_col(i, k, K)]) } } if (cs >= 4) { z = z \ TR[t] } if (cs >= 2) { z = z \ 1 } return(z) } // Regressor row for the generating y-model of unit i at time t. The // level terms (and the restricted deterministic term) are omitted // wherever the null has been imposed, so the column order matches the // restricted design ES[i].X exactly. // nodrift != 0 zeroes the deterministic entries of the row for an // equation whose null has been imposed. The deterministic coefficients // are still estimated in every replication, so their estimation cost is // reflected, but they are not propagated into the LEVEL of y*. // // This matters. Under the null the generating equation is // D.y* = c0 (+ c1 t) + sum_k gamma_k D.x* + u*, // so a fitted intercept makes y* a random walk WITH DRIFT and a fitted // trend makes it quadratic. The PSS (2001) Table CI null DGP, and the // authors' own OxGauss code, generate y* as a driftless random walk. // Propagating the fitted deterministics instead produces a null // distribution with a very heavy right tail for strongly trending // series -- critical values in the dozens rather than single digits -- // which is not the distribution the bounds test is calibrated against. real colvector _surardl_yrow(real matrix LVs, real colvector TR, struct surardl_eq vector E, real scalar i, real scalar K, real scalar t, real scalar cs, real scalar restricted, real scalar nodrift) { real scalar k, j, c0, ck real colvector z c0 = _surardl_col(i, 0, K) z = J(0, 1, .) if (restricted == 0) { z = z \ LVs[t-1, c0] for (k = 1; k <= K; k++) { z = z \ LVs[t-1, _surardl_col(i, k, K)] } } for (j = 1; j <= E[i].p - 1; j++) { z = z \ (LVs[t-j, c0] - LVs[t-j-1, c0]) } for (k = 1; k <= K; k++) { ck = _surardl_col(i, k, K) for (j = 0; j <= E[i].q[k] - 1; j++) { z = z \ (LVs[t-j, ck] - LVs[t-j-1, ck]) } } if (cs >= 4) { if (!(cs == 4 & restricted != 0)) { if (restricted != 0 & nodrift != 0) { z = z \ 0 } else { z = z \ TR[t] } } } if (cs >= 2) { if (!(cs == 2 & restricted != 0)) { if (restricted != 0 & nodrift != 0) { z = z \ 0 } else { z = z \ 1 } } } return(z) } // ================================================================== // Residual bootstrap for the SUR bounds statistics. // // Design: McNown, Sam & Goh (2018, section 3.2, steps 1-8), itself the // procedure of Palm, Smeekes & Urbain (2010, section 3.2 and Remark 2), // extended to the N-equation SUR system. This supplies the "Approach B" // that section 2.4 of Arif & Furuoka (2026) calls for but leaves // unspecified. // // Two features make the critical values SUR-specific rather than // single-equation: // * the N y-equation innovations and the N*K marginal innovations // are resampled as ONE time-indexed vector, so the contemporaneous // cross-equation covariance Sigma and the within-unit y-x // dependence are both preserved; // * every replication re-runs the ENTIRE statistic -- first-step // OLS, re-estimation of Sigma-hat, FGLS solve and Wald form -- so // the bootstrap distribution carries the same estimation noise as // the observed statistic. // // target == 0 : impose H0 on every y-equation at once (the null DGP of // the authors' own Monte Carlo code); one pass yields // critical values for all N equations and the joint test // target == i : impose H0 on equation i only, leaving the other y // equations and all marginal equations unrestricted // (faithful to McNown et al.; the default) // // dgp == 1 : marginal VECM of order mo fitted to the data // dgp == 0 : regressors generated as pure random walks, the DGP behind // the PSS and Narayan tabulations (default) // // nodrift != 0 : generate y* under the null as a DRIFTLESS random walk, // as in the PSS Table CI design. The default 0 propagates the // fitted intercept and trend into the level of y*, which is // the literal reading of McNown et al. step 4 and keeps the // bootstrap faithful to the data's own deterministic // behaviour. // // Returns B x 3N: columns 1..N are F*, N+1..2N are t*, 2N+1..3N are // Fidv*. // ================================================================== real matrix _surardl_boot(real matrix LV, real colvector TR, real scalar N, real scalar K, real matrix LAGS, real scalar cs, real scalar dfk, real scalar iterate, real scalar tol, real scalar B, real scalar target, real scalar dgp, real scalar mo, real scalar nodrift, string scalar dep, string rowvector xn, real scalar nfail) { struct surardl_eq vector E, EG, ES real scalar i, k, b, t, L, Lg, T, Tb, Tf, M, c0, ck, jr, rest, rc real scalar nit, conv real matrix Out, S, U0, V, URes, UR, VR, VRall, VRs, LVs, Bmarg, St real colvector ths, bols, seols, rssv, r2v, rmsev, keep, idx, zrow real rowvector np, dxs Tf = rows(LV) L = _surardl_lagmax(LAGS, K) T = Tf - L Lg = max((L, mo + 1)) Tb = Tf - Lg Out = J(B, 3 * N, .) nfail = 0 // ------- Step 1: generating models, null imposed where required -- E = _surardl_build(LV, TR, N, K, LAGS, cs, dep, xn) EG = E for (i = 1; i <= N; i++) { if (target == 0 | target == i) { keep = J(0, 1, .) for (k = 1; k <= E[i].m; k++) { if (!anyof(E[i].rst, k)) { keep = keep \ k } } EG[i].X = E[i].X[., keep] EG[i].m = rows(keep) } } M = 0 for (i = 1; i <= N; i++) { EG[i].off = M M = M + EG[i].m } U0 = _surardl_ols(EG, bols, seols, rssv, r2v, rmsev) _surardl_fgls(EG, U0, dfk, iterate, tol, ths, V, S, URes, nit, conv) // ------- marginal models for the regressors --------------------- Bmarg = J(0, 0, .) VRall = J(Tb, N * K, 0) if (dgp != 0) { for (i = 1; i <= N; i++) { VR = J(0, 0, .) St = _surardl_marg(LV, TR, i, K, Lg, mo, cs, VR) if (i == 1) { Bmarg = J(rows(St), N * K, 0) } Bmarg[., ((i-1)*K+1)..(i*K)] = St VRall[., ((i-1)*K+1)..(i*K)] = VR } } else { for (i = 1; i <= N; i++) { for (k = 1; k <= K; k++) { ck = _surardl_col(i, k, K) VRall[., (i-1)*K+k] = LV[((Lg+1)::Tf), ck] - LV[(Lg::(Tf-1)), ck] } } } // ------- Step 2: rescale and recentre, aligned on Tb rows ------- np = J(1, N, 0) for (i = 1; i <= N; i++) { np[i] = EG[i].m } UR = _surardl_rescale(URes[((T - Tb + 1)::T), .], np) VRs = _surardl_rescale(VRall, J(1, N * K, cols(Bmarg))) // ------- Steps 3-8 ---------------------------------------------- for (b = 1; b <= B; b++) { // Step 3: ONE index draw per period, shared by every equation. idx = ceil(runiform(Tb, 1) :* Tb) idx = idx :+ (idx :< 1) // Steps 4-5: generate pseudo-levels recursively. Rows up to Lg // are kept at their observed values as initial conditions. LVs = LV for (t = Lg + 1; t <= Tf; t++) { jr = idx[t - Lg] // regressors first: they do not depend on the current dy for (i = 1; i <= N; i++) { dxs = J(1, K, 0) if (dgp != 0) { zrow = _surardl_margrow(LVs, TR, i, K, t, mo, cs) } for (k = 1; k <= K; k++) { ck = _surardl_col(i, k, K) if (dgp == 0) { dxs[k] = VRs[jr, (i-1)*K+k] } else { dxs[k] = (zrow' * Bmarg[., (i-1)*K+k]) + VRs[jr, (i-1)*K+k] } LVs[t, ck] = LVs[t-1, ck] + dxs[k] } } // dependent variables, using the just-generated dx for (i = 1; i <= N; i++) { c0 = _surardl_col(i, 0, K) rest = 0 if (target == 0 | target == i) { rest = 1 } zrow = _surardl_yrow(LVs, TR, E, i, K, t, cs, rest, nodrift) LVs[t, c0] = LVs[t-1, c0] + (zrow' * ths[|EG[i].off+1 \ EG[i].off+EG[i].m|]) + UR[jr, i] } } // Explosion guard. An unrestricted marginal VECM estimated on a // short sample can have roots outside the unit circle, in which // case the recursion diverges and the replication is worthless. // Such draws are discarded and counted, never used. rc = 0 if (hasmissing(LVs)) { rc = 1 } else { if (max(abs(LVs)) > 1e4 * (1 + max(abs(LV)))) { rc = 1 } } if (rc != 0) { nfail = nfail + 1 continue } // Steps 6-7: re-run the whole statistic on the pseudo-data. ES = _surardl_build(LVs, TR, N, K, LAGS, cs, dep, xn) St = J(N, 3, .) rc = _surardl_try(ES, dfk, iterate, tol, St) if (rc == 0) { Out[b, (1..N)] = St[., 1]' Out[b, ((N+1)..(2*N))] = St[., 2]' Out[b, ((2*N+1)..(3*N))] = St[., 3]' } else { nfail = nfail + 1 } } return(Out) } // Empirical quantile of a column (linear interpolation on the order // statistics), used for bootstrap and simulated critical values. real scalar _surardl_quant(real colvector x, real scalar p) { real colvector s real scalar n, h s = sort(select(x, x :< .), 1) n = rows(s) if (n < 2) { return(.) } h = p * n if (h < 1) { return(s[1]) } if (h >= n) { return(s[n]) } return(s[floor(h)] + (h - floor(h)) * (s[floor(h)+1] - s[floor(h)])) } // Bootstrap p-value (1 + #{beyond}) / (B + 1). // upper != 0 for right-tailed statistics, upper == 0 for the t statistic. real scalar _surardl_bootp(real colvector x, real scalar stat, real scalar upper) { real colvector s real scalar n s = select(x, x :< .) n = rows(s) if (n < 1) { return(.) } if (upper != 0) { return((1 + colsum(s :>= stat)) / (n + 1)) } return((1 + colsum(s :<= stat)) / (n + 1)) } // ================================================================== // Approach A of section 2.4: SUR-specific Monte Carlo critical value // BOUNDS. Data are generated under H0 with the estimated Sigma-hat as // the y-innovation covariance and the regressors either purely I(1) // (upper bound) or purely I(0) (lower bound), as in the notes to PSS // (2001, Table CI) and in the authors' own OxGauss // empirical-distribution code. // // i1 != 0 : x_t = x_{t-1} + e_t (upper bound) // i1 == 0 : x_t = e_t (lower bound) // rhox : equicorrelation imposed across units on the regressor // innovations. rhox = 0 reproduces the PSS design (the // regressors are independent across units); rhox = 1 // reproduces the authors' OxGauss design, in which every // unit's regressors are driven by the same shocks. // Theorem 3(c) makes the null distribution depend on the // other equations' regressors, so this is a real knob. // // Returns B x 3N, laid out as in _surardl_boot(). // ================================================================== real matrix _surardl_sim(real matrix S, real scalar N, real scalar K, real matrix LAGS, real scalar cs, real scalar dfk, real scalar iterate, real scalar tol, real scalar B, real scalar T, real scalar i1, real scalar burn, real scalar rhox, string scalar dep, string rowvector xn, real scalar nfail) { struct surardl_eq vector ES real scalar b, i, k, t, L, Tf, Tt, rc, sr real matrix Out, C, EPS, LVs, St, EX, CM real colvector TR L = _surardl_lagmax(LAGS, K) Tf = T + L Tt = Tf + burn Out = J(B, 3 * N, .) nfail = 0 C = cholesky(S)' TR = (1 :: Tf) // equicorrelation factor for the regressor innovations sr = sqrt(abs(rhox)) if (rhox < 0) { sr = 0 } for (b = 1; b <= B; b++) { EPS = invnormal(runiform(Tt, N)) * C EX = invnormal(runiform(Tt, N * K)) if (sr > 0) { CM = invnormal(runiform(Tt, K)) for (i = 1; i <= N; i++) { for (k = 1; k <= K; k++) { EX[., (i-1)*K+k] = sr :* CM[., k] + sqrt(1 - sr * sr) :* EX[., (i-1)*K+k] } } } LVs = J(Tt, N * (K + 1), 0) for (t = 2; t <= Tt; t++) { for (i = 1; i <= N; i++) { LVs[t, _surardl_col(i, 0, K)] = LVs[t-1, _surardl_col(i, 0, K)] + EPS[t, i] for (k = 1; k <= K; k++) { if (i1 != 0) { LVs[t, _surardl_col(i, k, K)] = LVs[t-1, _surardl_col(i, k, K)] + EX[t, (i-1)*K+k] } else { LVs[t, _surardl_col(i, k, K)] = EX[t, (i-1)*K+k] } } } } LVs = LVs[((burn+1)::Tt), .] ES = _surardl_build(LVs, TR, N, K, LAGS, cs, dep, xn) St = J(N, 3, .) rc = _surardl_try(ES, dfk, iterate, tol, St) if (rc == 0) { Out[b, (1..N)] = St[., 1]' Out[b, ((N+1)..(2*N))] = St[., 2]' Out[b, ((2*N+1)..(3*N))] = St[., 3]' } else { nfail = nfail + 1 } } return(Out) } // ================================================================ // from _sa_tests.mata // ================================================================ // ================================================================== // Joint / subset bounds restriction matrix, equations (65)-(67). // units is a row vector of equation indices to include; passing // (1..N) gives R_Joint of equation (67). // ================================================================== real matrix _surardl_Rjoint(struct surardl_eq vector E, real rowvector units, real scalar M) { real scalar u real matrix R R = J(0, M, .) for (u = 1; u <= length(units); u++) { R = R \ _surardl_R(E, units[u], M) } return(R) } // ================================================================== // Poolability restriction matrix, section 4.5(a): // H0: beta_{k,i} = beta_{k,1}, i = 2..N, k = 1..K // with K(N-1) restrictions. what selects the parameter block: // 1 = long-run level coefficients beta (K(N-1) restrictions) // 2 = speed of adjustment lambda (N-1 restrictions) // 3 = short-run coefficients gamma/psi (requires equal lag orders) // 4 = beta and lambda jointly // Returns J(0, M, .) when the restrictions are not well defined. // ================================================================== real matrix _surardl_Rpool(struct surardl_eq vector E, real scalar M, real scalar K, real scalar what) { real scalar i, k, N, ns real matrix R real rowvector r N = length(E) R = J(0, M, .) if (N < 2) { return(R) } if (what == 1 | what == 4) { for (i = 2; i <= N; i++) { for (k = 1; k <= K; k++) { r = J(1, M, 0) r[E[i].off + E[i].ibeta[k]] = 1 r[E[1].off + E[1].ibeta[k]] = -1 R = R \ r } } } if (what == 2 | what == 4) { for (i = 2; i <= N; i++) { r = J(1, M, 0) r[E[i].off + E[i].ilam] = 1 r[E[1].off + E[1].ilam] = -1 R = R \ r } } if (what == 3) { ns = length(E[1].ipsi) + length(E[1].igam) for (i = 2; i <= N; i++) { if (length(E[i].ipsi) + length(E[i].igam) != ns) { return(J(0, M, .)) } } for (i = 2; i <= N; i++) { for (k = 1; k <= length(E[1].ipsi); k++) { r = J(1, M, 0) r[E[i].off + E[i].ipsi[k]] = 1 r[E[1].off + E[1].ipsi[k]] = -1 R = R \ r } for (k = 1; k <= length(E[1].igam); k++) { r = J(1, M, 0) r[E[i].off + E[i].igam[k]] = 1 r[E[1].off + E[1].igam[k]] = -1 R = R \ r } } } return(R) } // Pairwise restriction on the speed of adjustment, section 4.5(c): // H0: lambda_i = lambda_j real matrix _surardl_Rlam(struct surardl_eq vector E, real scalar i, real scalar j, real scalar M) { real rowvector r r = J(1, M, 0) r[E[i].off + E[i].ilam] = 1 r[E[j].off + E[j].ilam] = -1 return(r) } // ================================================================== // Gradient of the long-run coefficients with respect to the system // parameter vector. Row (i-1)*K+k holds d theta_{k,i} / d theta, with // d theta_{k,i} / d lambda_i = beta_{k,i} / lambda_i^2 // d theta_{k,i} / d beta_{k,i} = -1 / lambda_i // This is the multivariate delta method of Remark 10, equations // (92)-(93), written once for the whole system so that ANY linear // hypothesis C * theta^LR = c can be tested as // W = (C th - c)' [C G V G' C']^-1 (C th - c). // ================================================================== real matrix _surardl_lrgrad(struct surardl_eq vector E, real colvector th, real scalar K, real scalar M) { real scalar i, k, N, il, ib, lam, bet real matrix G N = length(E) G = J(N * K, M, 0) for (i = 1; i <= N; i++) { il = E[i].off + E[i].ilam lam = th[il] for (k = 1; k <= K; k++) { ib = E[i].off + E[i].ibeta[k] bet = th[ib] G[(i-1)*K + k, il] = bet / (lam * lam) G[(i-1)*K + k, ib] = -1 / lam } } return(G) } real colvector _surardl_lrvec(struct surardl_eq vector E, real colvector th, real scalar K) { real scalar i, k, N real colvector v N = length(E) v = J(N * K, 1, .) for (i = 1; i <= N; i++) { for (k = 1; k <= K; k++) { v[(i-1)*K + k] = -th[E[i].off + E[i].ibeta[k]] / th[E[i].off + E[i].ilam] } } return(v) } // Wald test of C * theta^LR = 0. Returns (W, df, p). real rowvector _surardl_lrwald(real matrix C, real colvector thLR, real matrix G, real matrix V) { real colvector d real matrix Q real scalar w, df d = C * thLR Q = C * G * V * G' * C' w = (d' * invsym(Q) * d) df = rank(C) return((w, df, chi2tail(df, w))) } // Contrast selecting theta_{k,i} - theta_{k,j} for one regressor k. real matrix _surardl_Clrpair(real scalar i, real scalar j, real scalar k, real scalar N, real scalar K) { real rowvector c c = J(1, N * K, 0) c[(i-1)*K + k] = 1 c[(j-1)*K + k] = -1 return(c) } // Contrast selecting all K long-run differences between units i and j. real matrix _surardl_Clrpairall(real scalar i, real scalar j, real scalar N, real scalar K) { real scalar k real matrix C C = J(0, N * K, .) for (k = 1; k <= K; k++) { C = C \ _surardl_Clrpair(i, j, k, N, K) } return(C) } // Contrast for H0: theta_{k,1} = ... = theta_{k,N} (N-1 restrictions). real matrix _surardl_Clrpool(real scalar k, real scalar N, real scalar K) { real scalar i real matrix C C = J(0, N * K, .) for (i = 2; i <= N; i++) { C = C \ _surardl_Clrpair(i, 1, k, N, K) } return(C) } // ================================================================== // Multiple-testing adjustment of a vector of p-values. // method 1 = Bonferroni, 2 = Holm, 3 = Sidak (all capped at 1). // Section 4.5 of the paper recommends such a correction for the // pairwise long-run equality tests, where the maximum pairwise // rejection rate reaches 53% at N = 5, T = 50. // ================================================================== real colvector _surardl_padj(real colvector p, real scalar method) { real scalar m, i, mx real colvector a, o, s m = rows(p) a = p if (m < 1) { return(a) } if (method == 1) { a = p :* m } if (method == 3) { a = 1 :- (1 :- p) :^ m } if (method == 2) { o = order(p, 1) s = p[o] mx = 0 for (i = 1; i <= m; i++) { s[i] = s[i] * (m - i + 1) if (s[i] < mx) { s[i] = mx } mx = s[i] } a[o] = s } for (i = 1; i <= m; i++) { if (a[i] > 1) { a[i] = 1 } } return(a) } // ================================================================== // Sigma-hat quality diagnostics. Returns // (N, T, T/N, reciprocal condition number, min eigenvalue, // 1 if positive definite else 0, mean |off-diagonal correlation|, // relative Frobenius norm of the off-diagonal part) // The last quantity is the natural reading of the "RFN" referred to in // section 5 of the paper; its defining equation is not given in the // draft, so it is documented here as surardl's own definition: // RFN = ||Sigma - diag(Sigma)||_F / ||Sigma||_F. // ================================================================== real rowvector _surardl_sigdiag(real matrix S, real scalar T) { real scalar N, i, j, pd, sa, nf, df, fo, ft real matrix R real colvector ev N = rows(S) R = S :/ sqrt(diagonal(S) * diagonal(S)') ev = symeigenvalues(S)' pd = 1 if (min(ev) <= 0) { pd = 0 } sa = 0 df = N * (N - 1) / 2 for (i = 2; i <= N; i++) { for (j = 1; j <= i - 1; j++) { sa = sa + abs(R[i, j]) } } fo = 0 ft = 0 for (i = 1; i <= N; i++) { for (j = 1; j <= N; j++) { ft = ft + S[i, j] * S[i, j] if (i != j) { fo = fo + S[i, j] * S[i, j] } } } nf = sqrt(fo) / sqrt(ft) return((N, T, T / N, min(ev) / max(ev), min(ev), pd, sa / df, nf)) } // ================================================================ // from _sa_load.mata // ================================================================ // ================================================================== // Assemble the matrix of level series required by _surardl_build(). // // Reads the long panel (one row per unit-period) and lays it out as // LV (Tfull x N*(K+1)), with the column block of unit i ordered // y_i, x_1i, ..., x_Ki. The window [t0, t1] must already be balanced // and gap-free; the caller guarantees this. // // ids holds the panel values in the order the equations are to appear. // Returns 0 on success, or the 1-based row of the first cell that // could not be filled (which would indicate an unbalanced window). // ================================================================== real scalar _surardl_load(string scalar vlist, string scalar touse, string scalar ivar, string scalar tvar, real colvector ids, real scalar t0, real scalar delta, real scalar Tfull, real scalar K, string scalar lvname) { real matrix D, LV, FL real colvector uv, tv real scalar n, r, i, ti, v, N, c N = rows(ids) D = st_data(., vlist, touse) uv = st_data(., ivar, touse) tv = st_data(., tvar, touse) n = rows(D) LV = J(Tfull, N * (K + 1), .) FL = J(Tfull, N, 0) for (r = 1; r <= n; r++) { ti = (tv[r] - t0) / delta + 1 if (ti < 1 | ti > Tfull) { continue } i = 0 for (c = 1; c <= N; c++) { if (ids[c] == uv[r]) { i = c } } if (i == 0) { continue } for (v = 0; v <= K; v++) { LV[ti, _surardl_col(i, v, K)] = D[r, v + 1] } FL[ti, i] = 1 } if (min(FL) == 0) { return(1) } st_matrix(lvname, LV) return(0) } // Lay a user-supplied trend variable out on the estimation window. void _sa_trendvar(string scalar v, string scalar touse, string scalar tvar, real scalar t0, real scalar delta, real scalar Tfull, string scalar trname) { real colvector tv, xv, TR real scalar r, ti xv = st_data(., v, touse) tv = st_data(., tvar, touse) TR = J(Tfull, 1, .) for (r = 1; r <= rows(tv); r++) { ti = (tv[r] - t0) / delta + 1 if (ti >= 1 & ti <= Tfull) { TR[ti] = xv[r] } } st_matrix(trname, TR) } // ================================================================ // from _sa_main.mata // ================================================================ // ================================================================== // Information-criterion lag selection, Remark 2 and equations // (19)-(20). Every candidate is fitted on the SAME effective sample, // the one implied by the maximum lag orders, so the criteria are // comparable across candidates and Sigma-hat is defined on a common // window. AIC = T ln(s2) + 2m and SBC = T ln(s2) + m ln T; these are // monotone equivalents of the paper's equations (19)-(20). // crit: 1 = AIC, 2 = SBC. // Returns the N x (1+K) matrix of selected orders. // ================================================================== // Helper: unit i keeps its candidate orders while every other unit is // pinned at the maximum orders, so _surardl_lagmax() -- and therefore // the effective sample -- stays fixed across candidates. real matrix _surardl_padlags(real matrix LAGS, real matrix LFIX, real scalar i, real scalar K) { real matrix L L = LFIX L[i, .] = LAGS[i, .] return(L) } real matrix _surardl_lagsel(real matrix LV, real colvector TR, real scalar N, real scalar K, real scalar cs, real scalar pmax, real rowvector qmax, real scalar crit, real scalar qmin, real scalar maxcomb, real scalar ncomb) { real matrix LAGS, LFIX, cand real scalar i, j, k, p, nc, T, bic, s2, m, v, ok real rowvector q, qb struct surardl_eq vector E real colvector bols, seols, rssv, r2v, rmsev real matrix U // enumerate candidate (p, q_1..q_K) combinations nc = pmax for (k = 1; k <= K; k++) { nc = nc * (qmax[k] - qmin + 1) } ncomb = nc LAGS = J(N, K + 1, 1) if (nc > maxcomb) { return(J(0, 0, .)) } cand = J(nc, K + 1, .) j = 0 for (p = 1; p <= pmax; p++) { q = J(1, K, qmin) for (v = 1; v <= nc / pmax; v++) { j = j + 1 cand[j, 1] = p cand[j, (2..(K+1))] = q // odometer increment over q k = K ok = 1 while (k >= 1 & ok) { q[k] = q[k] + 1 if (q[k] > qmax[k]) { q[k] = qmin k = k - 1 } else { ok = 0 } } } } // the fixed sample is set by the maximum orders LFIX = J(N, K + 1, 1) LFIX[., 1] = J(N, 1, pmax) for (k = 1; k <= K; k++) { LFIX[., k+1] = J(N, 1, qmax[k]) } T = rows(LV) - _surardl_lagmax(LFIX, K) for (i = 1; i <= N; i++) { bic = . qb = LAGS[i, .] for (j = 1; j <= nc; j++) { LAGS[i, .] = cand[j, .] // build only unit i on the fixed window E = _surardl_build(LV, TR, N, K, _surardl_padlags(LAGS, LFIX, i, K), cs, "y", J(1, K, "x")) if (rows(E[i].y) != T) { continue } if (rank(quadcross(E[i].X, E[i].X)) < E[i].m) { continue } U = _surardl_ols(E, bols, seols, rssv, r2v, rmsev) m = E[i].m s2 = rssv[i] / T if (crit == 1) { v = T * ln(s2) + 2 * m } else { v = T * ln(s2) + m * ln(T) } if (v < bic) { bic = v qb = cand[j, .] } } LAGS[i, .] = qb } return(LAGS) } // ================================================================== // Main driver. Reads its configuration from Stata, estimates the // system, runs the requested inference and posts everything back. // ================================================================== void _surardl_main() { struct surardl_eq vector E real matrix LV, LAGS, S, Si, U, U0, V, G, R, C, B1, B0 real matrix bnd, cvF, cvt, cvN, decF, decN, dect, LR, FI, EFF, FIT real matrix bcv, bpv, scv0, scv1, SIMd string matrix stripe real colvector TR, th, bols, seols, rssv, r2v, rmsev, thLR, Fr real rowvector w, dg, bp, sd, qmax real scalar N, K, cs, dfk, iterate, tol, lev, M, T, i, k, j, z, il, ib real scalar nit, conv, breps, bseed, btarget, bdgp, blags, sreps, srho, bndrift real scalar bglag, rssform, nfail, pmax, crit, maxcomb, ncomb, qmin real scalar lam, bet, vll, vbb, vlb, nf2 string scalar dep string rowvector xn string matrix eqn LV = st_matrix("__sa_LV") TR = st_matrix("__sa_TR") N = st_numscalar("__sa_N") K = st_numscalar("__sa_K") cs = st_numscalar("__sa_case") dfk = st_numscalar("__sa_dfk") iterate = st_numscalar("__sa_iterate") tol = st_numscalar("__sa_tol") lev = st_numscalar("__sa_level") rssform = st_numscalar("__sa_rssform") bglag = st_numscalar("__sa_bglags") breps = st_numscalar("__sa_breps") bseed = st_numscalar("__sa_bseed") btarget = st_numscalar("__sa_btarget") bdgp = st_numscalar("__sa_bdgp") blags = st_numscalar("__sa_blags") bndrift = st_numscalar("__sa_bndrift") sreps = st_numscalar("__sa_sreps") srho = st_numscalar("__sa_srho") dep = st_local("depvar") xn = tokens(st_local("xvars")) eqn = tokens(st_local("eqnames"))' // ---------------- lag orders ---------------- if (st_numscalar("__sa_dosel") != 0) { pmax = st_numscalar("__sa_pmax") qmax = st_matrix("__sa_qmax") crit = st_numscalar("__sa_crit") qmin = st_numscalar("__sa_qmin") maxcomb = st_numscalar("__sa_maxcomb") ncomb = 0 LAGS = _surardl_lagsel(LV, TR, N, K, cs, pmax, qmax, crit, qmin, maxcomb, ncomb) st_numscalar("__sa_ncomb", ncomb) if (rows(LAGS) == 0) { st_numscalar("__sa_rc", 459) return } } else { LAGS = st_matrix("__sa_lags") } st_matrix("__sa_lagsout", LAGS) // ---------------- estimation ---------------- E = _surardl_build(LV, TR, N, K, LAGS, cs, dep, xn) M = _surardl_M(E) T = rows(E[1].y) st_numscalar("__sa_T", T) st_numscalar("__sa_M", M) for (i = 1; i <= N; i++) { if (T - E[i].m < 1) { st_numscalar("__sa_rc", 2001) return } if (rank(quadcross(E[i].X, E[i].X)) < E[i].m) { st_numscalar("__sa_rc", 459) st_numscalar("__sa_badeq", i) return } } if (T <= N) { st_numscalar("__sa_rc", 2001) return } U0 = _surardl_ols(E, bols, seols, rssv, r2v, rmsev) if (st_numscalar("__sa_ols") != 0) { // equation-by-equation OLS benchmark: block-diagonal variance S = _surardl_sigma(U0, E, dfk) th = bols V = J(M, M, 0) for (i = 1; i <= N; i++) { V[|E[i].off+1, E[i].off+1 \ E[i].off+E[i].m, E[i].off+E[i].m|] = (rssv[i] / (T - E[i].m)) * invsym(quadcross(E[i].X, E[i].X)) } U = U0 nit = 0 conv = 1 } else { _surardl_fgls(E, U0, dfk, iterate, tol, th, V, S, U, nit, conv) } Si = invsym(S) st_numscalar("__sa_nit", nit) st_numscalar("__sa_conv", conv) st_matrix("__sa_Sigma", S) st_matrix("__sa_Corr", S :/ sqrt(diagonal(S) * diagonal(S)')) // coefficient stripe: equation name, coefficient name stripe = J(M, 2, "") for (i = 1; i <= N; i++) { for (k = 1; k <= E[i].m; k++) { stripe[E[i].off + k, 1] = eqn[i] stripe[E[i].off + k, 2] = E[i].cn[k] } } st_matrix("__sa_b", th') st_matrixcolstripe("__sa_b", stripe) st_matrix("__sa_V", V) st_matrixcolstripe("__sa_V", stripe) st_matrixrowstripe("__sa_V", stripe) // OLS standard errors and the efficiency gain of equation (95) EFF = J(M, 3, .) EFF[., 1] = seols EFF[., 2] = sqrt(diagonal(V)) EFF[., 3] = 1 :- EFF[., 2] :/ EFF[., 1] st_matrix("__sa_eff", EFF) st_matrixrowstripe("__sa_eff", stripe) // per-equation fit FIT = J(N, 5, .) for (i = 1; i <= N; i++) { FIT[i, 1] = quadcross(U[., i], U[., i]) FIT[i, 2] = rssv[i] FIT[i, 3] = r2v[i] FIT[i, 4] = T - E[i].m FIT[i, 5] = E[i].m } st_matrix("__sa_fit", FIT) // ---------------- bounds statistics ---------------- bnd = J(N, 5, .) for (i = 1; i <= N; i++) { R = _surardl_R(E, i, M) w = _surardl_wald(R, th, V) bnd[i, 1] = w[3] bnd[i, 2] = w[1] bnd[i, 4] = w[2] il = E[i].off + E[i].ilam bnd[i, 3] = th[il] / sqrt(V[il, il]) R = _surardl_Rbeta(E, i, M) w = _surardl_wald(R, th, V) bnd[i, 5] = w[3] } if (rssform != 0) { Fr = _surardl_rssF(E, Si, FIT[., 1], rssform - 1) st_matrix("__sa_Frss", Fr) } st_matrix("__sa_bnd", bnd) // joint test, equations (78)-(79) R = _surardl_Rjoint(E, (1..N), M) w = _surardl_wald(R, th, V) st_matrix("__sa_joint", (w[1], w[2], w[3], chi2tail(w[2], w[1]))) // tabulated critical values (k and the case are common to all // equations, so one lookup serves the whole system) cvF = _surardl_cvF_pss(cs, K) cvt = _surardl_cvt_pss(cs, K) cvN = _surardl_cvF_nar(cs, K, T) st_matrix("__sa_cvF", cvF) st_matrix("__sa_cvt", cvt) st_matrix("__sa_cvN", cvN) decF = J(N, 4, .) dect = J(N, 4, .) decN = J(N, 3, .) for (i = 1; i <= N; i++) { for (j = 1; j <= 4; j++) { decF[i, j] = _surardl_decide(bnd[i, 1], cvF[j, 1], cvF[j, 2], 1) dect[i, j] = _surardl_decide(bnd[i, 3], cvt[j, 1], cvt[j, 2], 0) } for (j = 1; j <= 3; j++) { decN[i, j] = _surardl_decide(bnd[i, 1], cvN[j, 1], cvN[j, 2], 1) } } st_matrix("__sa_decF", decF) st_matrix("__sa_dect", dect) st_matrix("__sa_decN", decN) // ---------------- long run ---------------- z = invnormal(1 - (1 - lev / 100) / 2) thLR = _surardl_lrvec(E, th, K) G = _surardl_lrgrad(E, th, K, M) LR = J(N * K, 8, .) FI = J(N * K, 3, .) for (i = 1; i <= N; i++) { il = E[i].off + E[i].ilam lam = th[il] vll = V[il, il] for (k = 1; k <= K; k++) { ib = E[i].off + E[i].ibeta[k] bet = th[ib] vbb = V[ib, ib] vlb = V[il, ib] j = (i-1)*K + k dg = G[j, .] LR[j, 1] = i LR[j, 2] = k LR[j, 3] = thLR[j] LR[j, 4] = sqrt(dg * V * dg') LR[j, 5] = LR[j, 3] / LR[j, 4] LR[j, 6] = 2 * normal(-abs(LR[j, 5])) LR[j, 7] = LR[j, 3] - z * LR[j, 4] LR[j, 8] = LR[j, 3] + z * LR[j, 4] FI[j, .] = _surardl_fieller(lam, bet, vll, vbb, vlb, z) } } st_matrix("__sa_lr", LR) st_matrix("__sa_fieller", FI) // ---------------- diagnostics ---------------- bp = _surardl_bpcd(U0) st_matrix("__sa_bpcd", bp) sd = _surardl_sigdiag(S, T) st_matrix("__sa_sigdiag", sd) if (bglag > 0) { st_matrix("__sa_bg", _surardl_bg(E, U, bglag)) } // residuals for predict and for postestimation reuse st_matrix("__sa_resid", U) st_matrix("__sa_residols", U0) // ---------------- bootstrap, Approach B ---------------- if (breps > 0) { rseed(bseed) if (btarget == 0) { B1 = _surardl_boot(LV, TR, N, K, LAGS, cs, dfk, iterate, tol, breps, 0, bdgp, blags, bndrift, dep, xn, nfail) } else { B1 = J(breps, 3 * N, .) nfail = 0 for (i = 1; i <= N; i++) { nf2 = 0 B0 = _surardl_boot(LV, TR, N, K, LAGS, cs, dfk, iterate, tol, breps, i, bdgp, blags, bndrift, dep, xn, nf2) B1[., i] = B0[., i] B1[., N + i] = B0[., N + i] B1[., 2 * N + i] = B0[., 2 * N + i] nfail = nfail + nf2 } } bcv = J(N, 9, .) bpv = J(N, 3, .) for (i = 1; i <= N; i++) { bcv[i, 1] = _surardl_quant(B1[., i], 0.90) bcv[i, 2] = _surardl_quant(B1[., i], 0.95) bcv[i, 3] = _surardl_quant(B1[., i], 0.99) bcv[i, 4] = _surardl_quant(B1[., N+i], 0.10) bcv[i, 5] = _surardl_quant(B1[., N+i], 0.05) bcv[i, 6] = _surardl_quant(B1[., N+i], 0.01) bcv[i, 7] = _surardl_quant(B1[., 2*N+i], 0.90) bcv[i, 8] = _surardl_quant(B1[., 2*N+i], 0.95) bcv[i, 9] = _surardl_quant(B1[., 2*N+i], 0.99) bpv[i, 1] = _surardl_bootp(B1[., i], bnd[i, 1], 1) bpv[i, 2] = _surardl_bootp(B1[., N+i], bnd[i, 3], 0) bpv[i, 3] = _surardl_bootp(B1[., 2*N+i], bnd[i, 5], 1) } st_matrix("__sa_bcv", bcv) st_matrix("__sa_bpv", bpv) st_numscalar("__sa_bfail", nfail) st_matrix("__sa_bdist", B1) } // ---------------- Monte Carlo bounds, Approach A ---------------- if (sreps > 0) { rseed(bseed) SIMd = _surardl_sim(S, N, K, LAGS, cs, dfk, iterate, tol, sreps, T, 1, 50, srho, dep, xn, nfail) scv1 = J(N, 3, .) for (i = 1; i <= N; i++) { scv1[i, 1] = _surardl_quant(SIMd[., i], 0.90) scv1[i, 2] = _surardl_quant(SIMd[., i], 0.95) scv1[i, 3] = _surardl_quant(SIMd[., i], 0.99) } nf2 = 0 SIMd = _surardl_sim(S, N, K, LAGS, cs, dfk, iterate, tol, sreps, T, 0, 50, srho, dep, xn, nf2) scv0 = J(N, 3, .) for (i = 1; i <= N; i++) { scv0[i, 1] = _surardl_quant(SIMd[., i], 0.90) scv0[i, 2] = _surardl_quant(SIMd[., i], 0.95) scv0[i, 3] = _surardl_quant(SIMd[., i], 0.99) } st_matrix("__sa_scv1", scv1) st_matrix("__sa_scv0", scv0) st_numscalar("__sa_sfail", nfail + nf2) } st_numscalar("__sa_rc", 0) } // ================================================================ // from _sa_estat.mata // ================================================================ // Rebuild the equation structures from what surardl left in e(). The // level series are stored in e(lv) so that every postestimation // command works on exactly the estimation window, whatever has // happened to the data in memory since. struct surardl_eq vector _surardl_reload() { real matrix LV, LAGS real colvector TR real scalar cs, N, K string scalar dep string rowvector xn LV = st_matrix("e(lv)") TR = st_matrix("e(tr)") LAGS = st_matrix("e(lags)") cs = st_numscalar("e(case)") N = rows(LAGS) K = cols(LAGS) - 1 dep = st_global("e(depvar)") xn = tokens(st_global("e(indepvars)")) return(_surardl_build(LV, TR, N, K, LAGS, cs, dep, xn)) } // ------------------------------------------------------------------ // Pairwise and joint equality of the speeds of adjustment, // H0: lambda_i = lambda_j and H0: lambda_1 = ... = lambda_N. // Section 4.5(c). Returns npairs x 5: i, j, chi2, p, adjusted p. // ------------------------------------------------------------------ void _sa_e_lamtest(real scalar padjm, string scalar outname, string scalar jointname) { struct surardl_eq vector E real matrix R, Out real colvector th, p real rowvector w real scalar N, M, i, j, r E = _surardl_reload() N = length(E) M = _surardl_M(E) th = st_matrix("e(b)")' Out = J(N * (N - 1) / 2, 5, .) r = 0 for (i = 1; i <= N - 1; i++) { for (j = i + 1; j <= N; j++) { r = r + 1 R = _surardl_Rlam(E, i, j, M) w = _surardl_wald(R, th, st_matrix("e(V)")) Out[r, 1] = i Out[r, 2] = j Out[r, 3] = w[1] Out[r, 4] = chi2tail(1, w[1]) } } p = Out[., 4] Out[., 5] = _surardl_padj(p, padjm) st_matrix(outname, Out) R = _surardl_Rpool(E, M, cols(st_matrix("e(lags)")) - 1, 2) w = _surardl_wald(R, th, st_matrix("e(V)")) st_matrix(jointname, (w[1], w[2], chi2tail(w[2], w[1]))) } // ------------------------------------------------------------------ // Pairwise and joint equality of the long-run coefficients for // regressor k, H0: theta_{k,i} = theta_{k,j}. Multivariate delta // method, equations (92)-(93) and Remark 10. // ------------------------------------------------------------------ void _sa_e_lrtest(real scalar k, real scalar padjm, string scalar outname, string scalar jointname) { struct surardl_eq vector E real matrix G, C, V, Out real colvector th, thLR, p real rowvector w real scalar N, M, K, i, j, r E = _surardl_reload() N = length(E) M = _surardl_M(E) K = cols(st_matrix("e(lags)")) - 1 th = st_matrix("e(b)")' V = st_matrix("e(V)") thLR = _surardl_lrvec(E, th, K) G = _surardl_lrgrad(E, th, K, M) Out = J(N * (N - 1) / 2, 6, .) r = 0 for (i = 1; i <= N - 1; i++) { for (j = i + 1; j <= N; j++) { r = r + 1 C = _surardl_Clrpair(i, j, k, N, K) w = _surardl_lrwald(C, thLR, G, V) Out[r, 1] = i Out[r, 2] = j Out[r, 3] = thLR[(i-1)*K+k] - thLR[(j-1)*K+k] Out[r, 4] = w[1] Out[r, 5] = w[3] } } p = Out[., 5] Out[., 6] = _surardl_padj(p, padjm) st_matrix(outname, Out) C = _surardl_Clrpool(k, N, K) w = _surardl_lrwald(C, thLR, G, V) st_matrix(jointname, w) } // ------------------------------------------------------------------ // Poolability tests, section 4.5(a). what: 1 = beta, 2 = lambda, // 3 = short-run, 4 = beta and lambda jointly. // ------------------------------------------------------------------ void _sa_e_pool(real scalar what, string scalar outname) { struct surardl_eq vector E real matrix R real colvector th real rowvector w real scalar M, K E = _surardl_reload() M = _surardl_M(E) K = cols(st_matrix("e(lags)")) - 1 th = st_matrix("e(b)")' R = _surardl_Rpool(E, M, K, what) if (rows(R) == 0) { st_matrix(outname, (., ., .)) return } w = _surardl_wald(R, th, st_matrix("e(V)")) st_matrix(outname, (w[1], w[2], chi2tail(w[2], w[1]))) } // Fieller sets recomputed at an arbitrary confidence level. void _sa_e_fieller(real scalar lev, string scalar outname) { struct surardl_eq vector E real matrix V, FI real colvector th real scalar N, K, i, k, il, ib, z E = _surardl_reload() N = length(E) K = cols(st_matrix("e(lags)")) - 1 th = st_matrix("e(b)")' V = st_matrix("e(V)") z = invnormal(1 - (1 - lev / 100) / 2) FI = J(N * K, 3, .) for (i = 1; i <= N; i++) { il = E[i].off + E[i].ilam for (k = 1; k <= K; k++) { ib = E[i].off + E[i].ibeta[k] FI[(i-1)*K+k, .] = _surardl_fieller(th[il], th[ib], V[il, il], V[ib, ib], V[il, ib], z) } } st_matrix(outname, FI) } // ------------------------------------------------------------------ // Percentile and bias-corrected bootstrap confidence intervals for // the long-run coefficients. The resampling is done under the // ALTERNATIVE (unrestricted residuals), which is the right null-free // design for an interval, and it preserves the cross-equation // covariance by drawing one time index for the whole system. // ------------------------------------------------------------------ void _sa_e_lrboot(real scalar B, real scalar seed, real scalar lev, string scalar outname) { struct surardl_eq vector E, ES real matrix LV, LAGS, V, U, U0, S, Out, D, LVs, Ub real colvector TR, th, bols, seols, rssv, r2v, rmsev, idx, zrow, ths real colvector thLR real rowvector np real scalar N, K, M, cs, T, Tf, L, i, k, t, b, c0, ck, nit, conv, dfk real scalar lo, hi, z0, nb string scalar dep string rowvector xn LV = st_matrix("e(lv)") TR = st_matrix("e(tr)") LAGS = st_matrix("e(lags)") cs = st_numscalar("e(case)") N = rows(LAGS) K = cols(LAGS) - 1 dep = st_global("e(depvar)") xn = tokens(st_global("e(indepvars)")) dfk = 0 E = _surardl_build(LV, TR, N, K, LAGS, cs, dep, xn) M = _surardl_M(E) Tf = rows(LV) L = _surardl_lagmax(LAGS, K) T = Tf - L U0 = _surardl_ols(E, bols, seols, rssv, r2v, rmsev) _surardl_fgls(E, U0, dfk, 0, 1e-9, ths, V, S, U, nit, conv) np = J(1, N, 0) for (i = 1; i <= N; i++) { np[i] = E[i].m } U = _surardl_rescale(U, np) D = J(B, N * K, .) rseed(seed) for (b = 1; b <= B; b++) { idx = ceil(runiform(T, 1) :* T) idx = idx :+ (idx :< 1) LVs = LV for (t = L + 1; t <= Tf; t++) { for (i = 1; i <= N; i++) { c0 = _surardl_col(i, 0, K) zrow = _surardl_yrow(LVs, TR, E, i, K, t, cs, 0, 0) LVs[t, c0] = LVs[t-1, c0] + (zrow' * ths[|E[i].off+1 \ E[i].off+E[i].m|]) + U[idx[t - L], i] } } ES = _surardl_build(LVs, TR, N, K, LAGS, cs, dep, xn) U0 = _surardl_ols(ES, bols, seols, rssv, r2v, rmsev) _surardl_fgls(ES, U0, dfk, 0, 1e-9, th, V, S, Ub, nit, conv) D[b, .] = _surardl_lrvec(ES, th, K)' } thLR = _surardl_lrvec(E, ths, K) lo = (1 - lev / 100) / 2 hi = 1 - lo Out = J(N * K, 5, .) for (k = 1; k <= N * K; k++) { Out[k, 1] = thLR[k] Out[k, 2] = _surardl_quant(D[., k], lo) Out[k, 3] = _surardl_quant(D[., k], hi) // bias-corrected percentile interval nb = colsum(select(D[., k], D[., k] :< .) :< thLR[k]) z0 = invnormal(nb / (B + 1)) Out[k, 4] = _surardl_quant(D[., k], normal(2 * z0 + invnormal(lo))) Out[k, 5] = _surardl_quant(D[., k], normal(2 * z0 + invnormal(hi))) } st_matrix(outname, Out) } end