*! bootur 1.0.0 20jul2026 *! Bootstrap Unit Root Tests -- Stata port of the R package bootUR (Smeekes & Wilms) *! Author: Merwan Roudane (merwanroudane920@gmail.com) *! GitHub: https://github.com/merwanroudane *------------------------------------------------------------------------------- * Dispatcher. *------------------------------------------------------------------------------- program define bootur, rclass version 14.0 gettoken sub rest : 0, parse(" ,") local sub = strlower("`sub'") if ("`sub'"=="_check") { _bucheck `rest' } else if ("`sub'"=="_pval") { _bupval `rest' } else if ("`sub'"=="adf") { _bu_adf `rest' } else if ("`sub'"=="ur") { _bu_boottest ur `rest' } else if ("`sub'"=="fdr") { _bu_boottest fdr `rest' } else if ("`sub'"=="sqt") { _bu_boottest sqt `rest' } else if ("`sub'"=="panel") { _bu_boottest panel `rest' } else if ("`sub'"=="union") { _bu_single union `rest' } else if ("`sub'"=="bootadf") { _bu_single bootadf `rest' } else if ("`sub'"=="order") { _bu_order `rest' } else if ("`sub'"=="diff") { _bu_diff `rest' } else if ("`sub'"=="plotmiss") { _bu_plotmiss `rest' } else if ("`sub'"=="plotorder") { _bu_plotorder `rest' } else { di as error "unknown bootur subcommand: `sub'" di as error "valid: adf ur union bootadf fdr sqt panel order diff plotmiss plotorder" exit 198 } return add end program define _bucheck, rclass * Hidden harness: compute the deterministic ADF engine for one series/spec. version 14.0 syntax varname(numeric) [if] [in], DET(int) IC(int) SCALE(int) TWO(int) /// PMIN(int) PMAX(int) [QD(int 0)] marksample touse, novarlist tempname g t lag mata: _bu_check_one("`varlist'", "`touse'", `det', `ic', `scale', `two', /// `pmin', `pmax', `qd', "`g'", "`t'", "`lag'") return scalar gamma = `g' return scalar tstat = `t' return scalar lag = `lag' return local icvec "`icvec'" return local shmin "`shmin'" return local umiss "`umiss'" return local ysmiss "`ysmiss'" end program define _bupval, rclass version 14.0 syntax , STAT(real) CASE(string) NOBS(real) tempname p mata: st_numscalar("`p'", bu_fpval(`stat', "`case'", `nobs')) return scalar pval = `p' end *------------------------------------------------------------------------------- * Shared driver for boot_ur / boot_fdr / boot_sqt / boot_panel *------------------------------------------------------------------------------- program define _bu_boottest, rclass version 14.0 gettoken mode 0 : 0 syntax varlist(numeric) [if] [in] [, BOOTstrap(string) B(integer 1999) /// BLOCKlength(integer -1) AR(real -1) MINlag(integer 0) MAXlag(integer -1) /// CRITerion(string) SCale(integer 1) DETerministics(string) DETRend(string) /// UNION(integer 1) LEVel(real 0.05) STEPs(numlist) SEED(integer -1) /// Dataname(string) noDOTs ] local B "`b'" * ---- resolve mode ---- local modenum = cond("`mode'"=="ur",1, cond("`mode'"=="fdr",2, /// cond("`mode'"=="sqt",3,4))) * ---- bootstrap method ---- if ("`bootstrap'"=="") local bootstrap "AWB" local bootstrap = strupper("`bootstrap'") local boot = 0 foreach k in MBB BWB DWB AWB SB SWB { local ++boot if ("`bootstrap'"=="`k'") local bootcode = `boot' } if ("`bootcode'"=="") { di as error "bootstrap() must be one of MBB BWB DWB AWB SB SWB" exit 198 } * ---- information criterion ---- if ("`criterion'"=="") local criterion "MAIC" local criterion = strupper("`criterion'") local ic = 0 foreach k in AIC BIC MAIC MBIC { local ++ic if ("`criterion'"=="`k'") local iccode = `ic' } if ("`iccode'"=="") { di as error "criterion() must be one of AIC BIC MAIC MBIC" exit 198 } * ---- sample ---- marksample touse, novarlist markout `touse', strok qui count if `touse' local T = r(N) local N : word count `varlist' if (`T' < 10) { di as error "too few observations" exit 2001 } if (`level'*(`B'+1) < 1) { di as error "B too low to perform the test at significance level `level'" exit 198 } * ---- union vs deterministic specification ---- if (`union') { if ("`deterministics'"!="") /// di as text "note: deterministics() ignored for the union test" if ("`detrend'"!="") /// di as text "note: detrend() ignored for the union test" local dc = 2 local detr = 1 } else { if ("`deterministics'"=="") local deterministics "intercept" local dc = cond("`deterministics'"=="none",0, /// cond("`deterministics'"=="intercept",1, /// cond("`deterministics'"=="trend",2,-1))) if (`dc'==-1) { di as error "deterministics() must be none, intercept or trend" exit 198 } if ("`detrend'"=="") local detrend "OLS" local detrend = strupper("`detrend'") local detr = cond("`detrend'"=="OLS",1,cond("`detrend'"=="QD",2,-1)) if (`detr'==-1) { di as error "detrend() must be OLS or QD" exit 198 } } * ---- block length / AR parameter ---- if (`blocklength' < 0) local blocklength = round(1.75*`T'^(1/3)) if (`ar' < 0) local ar = 0.01^(1/`blocklength') * ---- max lag (Schwert-type rule, small-sample corrected) ---- if (`maxlag' < 0) { local maxlag = floor(round(12*(`T'/100)^(1/4)) /// - 7*max(1-`T'/50,0)*(`T'/100)^(1/4)) if (`bootcode'==1) { local pib = 1 - (1 - (1/(`T'-`blocklength'))^(ceil(`T'/`blocklength')-1))^`B' if (`pib' > 0.01) local maxlag = min(`maxlag', `blocklength') } } if (`maxlag' < `minlag') local maxlag = `minlag' * ---- missing-value / unbalancedness handling ---- tempname jflag mata: st_numscalar("`jflag'", bu_check_missing("`varlist'", "`touse'")) if (`jflag' == 1) { di as error "internal missing values detected inside the sample" exit 416 } local joint = 1 tempname allbal mata: st_numscalar("`allbal'", bu_all_balanced("`varlist'", "`touse'")) if (`allbal'==0) { if (inlist(`bootcode',1,5)) { if (`modenum'==1 & `N'>1) { local joint = 0 di as text /// "note: missing values cause the resampling bootstrap to be run for each series individually" } else { di as error /// "resampling-based bootstraps MBB and SB cannot handle unbalanced series" exit 416 } } } if (inlist(`bootcode',5,6) & `modenum'!=1) /// di as text "note: SB and SWB are only recommended for the ur test; see help bootur" * ---- steps for the sequential quantile test ---- tempname PVEC if (`modenum'==3) { _bu_makesteps `N' "`steps'" matrix `PVEC' = r(pvec) } else { matrix `PVEC' = (0) } * ---- run the engine ---- if (`seed' >= 0) set seed `seed' local showdots = cond("`dots'"=="nodots",0,1) if (`showdots') di as text "Bootstrapping (B = `B') " _continue tempname OPT matrix `OPT' = (`bootcode', `B', `blocklength', `ar', `iccode', `scale', /// 0.1, `minlag', `maxlag', `union', `dc', `detr', `level', `joint', /// `modenum', `showdots') mata: bu_run("`varlist'", "`touse'", st_matrix("`OPT'"), "`PVEC'") if (`showdots') di as text " done" local D = __bu_D if ("`dataname'"=="") local dataname "`varlist'" local detrname = cond(`detr'==1,"OLS","QD") local detname = cond(`dc'==0,"none",cond(`dc'==1,"intercept","trend")) _bu_report, mode(`modenum') union(`union') n(`N') d(`D') level(`level') /// boot(`bootstrap') crit(`criterion') detname(`detname') detrname(`detrname') /// blocklength(`blocklength') breps(`B') maxlag(`maxlag') minlag(`minlag') /// dataname(`dataname') series(`varlist') return add end program define _bu_makesteps, rclass * Build the (increasing) p-vector for BSQT from user 'steps' (units or quantiles) args N steps if ("`steps'"=="") { numlist "0/`N'" local steps "`r(numlist)'" } * units (all non-negative integers) vs quantiles (values in [0,1]) local isunit = 1 foreach s of numlist `steps' { if (`s' != int(`s') | `s' < 0) local isunit = 0 } tempname pv if (`isunit') { local pl "`steps'" } else { * quantiles in [0,1] -> integer group boundaries local pl "" foreach s of numlist `steps' { local pl "`pl' `=round(`s'*`N')'" } } * sort, unique, clamp to <=N, prepend 0 / append N numlist "`pl'", sort local pl "`r(numlist)'" local clean "" local prev "." foreach s of numlist `pl' { if (`s' <= `N' & "`s'"!="`prev'") { local clean "`clean' `s'" local prev "`s'" } } local first : word 1 of `clean' if (`first' > 0) local clean "0 `clean'" local last : word `=wordcount("`clean'")' of `clean' if (`last' < `N') local clean "`clean' `N'" matrix `pv' = J(1, wordcount("`clean'"), 0) local i = 0 foreach s of numlist `clean' { local ++i matrix `pv'[1,`i'] = `s' } return matrix pvec = `pv' end *------------------------------------------------------------------------------- * Output formatter + stored results for the bootstrap tests *------------------------------------------------------------------------------- program define _bu_report, rclass version 14.0 syntax , mode(int) union(int) n(int) d(int) level(real) boot(string) /// crit(string) detname(string) detrname(string) blocklength(int) /// breps(int) maxlag(int) minlag(int) dataname(string) series(string) tempname est stat pval dstat dpar dlag dpval matrix `est' = __bu_est matrix `stat' = __bu_stat matrix `pval' = __bu_pval matrix `dstat'= __bu_dstat matrix `dpar' = __bu_dpar matrix `dlag' = __bu_dlag matrix `dpval'= __bu_dpval local unitxt = cond(`union',"union ","") if (`mode'==1) local mtxt "`boot' bootstrap `unitxt'test on each individual series" else if (`mode'==2) local mtxt "`boot' bootstrap `unitxt'test with false discovery rate control" else if (`mode'==3) local mtxt "`boot' bootstrap sequential quantile `unitxt'test" else local mtxt "Panel `boot' bootstrap group-mean `unitxt'test" di "" di as text "`mtxt'" di as text "Data: " as result "`dataname'" as text " (T-based block length = " /// as result "`blocklength'" as text ", B = " as result "`breps'" as text ")" if (`mode'==4) { di as text "H0: all series have a unit root" _col(48) /// "Ha: some series are stationary" } else { di as text "H0: series has a unit root" _col(48) "Ha: series is stationary" } if (!`union') di as text "Deterministics: " as result "`detname'" /// as text " Detrending: " as result "`detrname'" di as text "{hline 73}" * ---------- mode-specific body ---------- if (`mode'==4) { * panel group-mean local gm = __bu_gm local gmp = __bu_gmp di as text %-30s "Group-mean statistic" as result %14.4f `gm' di as text %-30s "Bootstrap p-value" as result %14.4f `gmp' di as text "{hline 73}" return scalar statistic = `gm' return scalar p_value = `gmp' return scalar N_series = `n' } else if (`mode'==1) { * individual tests di as text %-16s "series" %14s "est.root" %14s "statistic" %14s "p-value" %9s "reject" di as text "{hline 73}" local rej "" forvalues i = 1/`n' { local nm : word `i' of `series' local e = `est'[`i',1] local s = `stat'[`i',1] local p = `pval'[`i',1] local star = cond(`p'<`level',"*"," ") local rej "`rej' `=(`p'<`level')'" if (`union') { di as text %-16s abbrev("`nm'",16) %14s " ." /// as result %14.4f `s' %14.4f `p' %8s "`star'" } else { di as text %-16s abbrev("`nm'",16) as result %14.4f (`e'+1) /// %14.4f `s' %14.4f `p' as text %8s "`star'" } } di as text "{hline 73}" di as text "* rejects H0 at the `=`level'*100'% level" if (`n'==1) { return scalar statistic = `stat'[1,1] return scalar p_value = `pval'[1,1] return scalar estimate = `est'[1,1] } else { matrix rownames `stat' = `series' matrix rownames `pval' = `series' matrix rownames `est' = `series' return matrix statistic = `stat' return matrix p_value = `pval' return matrix estimate = `est' } } else { * fdr / sqt: sequence of tests + rejections tempname rej seq matrix `rej' = __bu_rej matrix `seq' = __bu_seq if (`mode'==2) { di as text "Step-down sequence (false discovery rate control):" di as text %-10s "series" %16s "statistic" %16s "crit. value" di as text "{hline 44}" local ns = rowsof(`seq') forvalues i = 1/`ns' { local idx = `seq'[`i',1] local nm : word `idx' of `series' di as text %-10s abbrev("`nm'",10) as result /// %16.4f `seq'[`i',2] %16.4f `seq'[`i',3] } } else { di as text "Sequential quantile test steps:" di as text %6s "step" %10s "H0 #I(0)" %10s "H1 #I(0)" %14s "statistic" %12s "p-value" di as text "{hline 54}" local ns = rowsof(`seq') forvalues i = 1/`ns' { di as text %6.0f `i' %10.0f `seq'[`i',1] %10.0f `seq'[`i',2] /// as result %14.4f `seq'[`i',4] %12.4f `seq'[`i',5] } } di as text "{hline 54}" * rejections summary local nrej = 0 forvalues i = 1/`n' { if (`rej'[`i',1]==1) local ++nrej } di as text "Series for which H0 (unit root) is rejected: " as result "`nrej'" as text " of `n'" forvalues i = 1/`n' { if (`rej'[`i',1]==1) { local nm : word `i' of `series' di as result " `nm'" } } matrix rownames `rej' = `series' return matrix rejections = `rej' return matrix sequence = `seq' } matrix rownames `dstat' = `series' matrix rownames `dpar' = `series' matrix rownames `dlag' = `series' matrix rownames `dpval' = `series' return matrix indiv_stat = `dstat' return matrix indiv_est = `dpar' return matrix indiv_lag = `dlag' return matrix indiv_pval = `dpval' return local bootstrap "`boot'" return local criterion "`crit'" return scalar level = `level' return scalar block_length = `blocklength' return scalar B = `breps' end *------------------------------------------------------------------------------- * Placeholder subcommands (implemented later in this file) *------------------------------------------------------------------------------- program define _bu_single, rclass version 14.0 gettoken which 0 : 0 if ("`which'"=="union") _bu_boottest ur `0' union(1) else _bu_boottest ur `0' union(0) return add end *------------------------------------------------------------------------------- * adf : standard augmented Dickey-Fuller test with MacKinnon asymptotic p-value *------------------------------------------------------------------------------- program define _bu_adf, rclass version 14.0 syntax varlist(numeric) [if] [in] [, DETerministics(string) MINlag(integer 0) /// MAXlag(integer -1) CRITerion(string) SCale(integer 1) ONEstep Dataname(string) ] local N : word count `varlist' if (`N' > 1) { di as error "adf handles a single time series; use {bf:bootur ur} for several" exit 198 } if ("`deterministics'"=="") local deterministics "intercept" local dc = cond("`deterministics'"=="none",0, /// cond("`deterministics'"=="intercept",1, /// cond("`deterministics'"=="trend",2,-1))) if (`dc'==-1) { di as error "deterministics() must be none, intercept or trend" exit 198 } if ("`criterion'"=="") local criterion "MAIC" local criterion = strupper("`criterion'") local ic = cond("`criterion'"=="AIC",1,cond("`criterion'"=="BIC",2, /// cond("`criterion'"=="MAIC",3,cond("`criterion'"=="MBIC",4,-1)))) if (`ic'==-1) { di as error "criterion() must be one of AIC BIC MAIC MBIC" exit 198 } local two = cond("`onestep'"=="",1,0) marksample touse, novarlist markout `touse', strok qui count if `touse' local T = r(N) if (`maxlag' < 0) { local maxlag = round(12*(`T'/100)^(1/4) - 7*max(1-`T'/50,0)*(`T'/100)^(1/4)) } if (`maxlag' < `minlag') local maxlag = `minlag' tempname g t lag tt mata: bu_adf_run("`varlist'", "`touse'", `dc', `ic', `scale', `two', /// `minlag', `maxlag', "`g'", "`t'", "`lag'", "`tt'") local gamma = `g' local tstat = `t' local lag = `lag' local TT = `tt' * MacKinnon asymptotic p-value local case = cond(`dc'==0,"nc",cond(`dc'==1,"c","ct")) local Nreg = max(`TT'-`maxlag'-1, 20) tempname p mata: st_numscalar("`p'", bu_fpval(`tstat', "`case'", `Nreg')) local pval = `p' if ("`dataname'"=="") local dataname "`varlist'" local sttxt = cond(`two',"Two-step","One-step") di "" di as text "`sttxt' augmented Dickey-Fuller test (with `deterministics') on a single series" di as text "Data: " as result "`dataname'" di as text "H0: series has a unit root" _col(48) "Ha: series is stationary" di as text "{hline 60}" di as text %-24s "estimate of largest root" as result %14.5f (`gamma'+1) di as text %-24s "ADF test statistic" as result %14.5f `tstat' di as text %-24s "MacKinnon p-value" as result %14.5f `pval' di as text %-24s "selected lag length" as result %14.0f `lag' di as text "{hline 60}" return scalar p_value = `pval' return scalar statistic= `tstat' return scalar estimate = `gamma' return scalar lag = `lag' return scalar N = `TT' return local criterion "`criterion'" return local deterministics "`deterministics'" end *------------------------------------------------------------------------------- * diff : difference multiple series, each by its own order (diff_mult) *------------------------------------------------------------------------------- program define _bu_diff, rclass version 14.0 syntax varlist(numeric) [if] [in], Orders(numlist integer >=0) /// [ GENerate(string) replace ] local N : word count `varlist' local no : word count `orders' if (`no' != `N') { di as error "orders() must have `N' values, one per variable" exit 198 } if ("`generate'"=="") local generate "d_" marksample touse, novarlist markout `touse', strok local outlist "" local i = 0 foreach v of local varlist { local ++i local d : word `i' of `orders' local nv "`generate'`v'" if ("`replace'"!="") capture drop `nv' qui gen double `nv' = . mata: bu_diff_store("`v'", "`nv'", "`touse'", `d') label variable `nv' "D`d'.`v'" local outlist "`outlist' `nv'" } di as text "created differenced series: " as result "`outlist'" return local newvars "`outlist'" end *------------------------------------------------------------------------------- * order : determine the order of integration of each series (order_integration) *------------------------------------------------------------------------------- program define _bu_order, rclass version 14.0 syntax varlist(numeric) [if] [in], [ Method(string) MAXOrder(integer 2) /// LEVel(real 0.05) GENerate(string) * ] if ("`method'"=="") local method "ur" local N : word count `varlist' marksample touse, novarlist markout `touse', strok tempname order matrix `order' = J(1, `N', .) * indices of series still under consideration local active "" forvalues i = 1/`N' { local active "`active' `i'" } forvalues di = `=`maxorder'-1'(-1)0 { local nact : word count `active' if (`nact'==0) continue, break * difference the active series di times into tempvars local tvlist "" local avars "" foreach idx of local active { local v : word `idx' of `varlist' local avars "`avars' `v'" tempvar tv qui gen double `tv' = . mata: bu_diff_store("`v'", "`tv'", "`touse'", `di') local tvlist "`tvlist' `tv'" } * run the chosen test and obtain a rejection vector quietly _bu_order_test `method' `tvlist' , level(`level') `options' matrix _rej = r(rej) * classify local newactive "" local k = 0 foreach idx of local active { local ++k local rj = _rej[1,`k'] if (`rj'==0) { matrix `order'[1,`idx'] = `di' + 1 } else { if (`di'==0) matrix `order'[1,`idx'] = 0 else local newactive "`newactive' `idx'" } } local active "`newactive'" } * display di "" di as text "Order of integration (method = `method', level = `level')" di as text "{hline 40}" di as text %-24s "series" %10s "order" di as text "{hline 40}" forvalues i = 1/`N' { local v : word `i' of `varlist' di as text %-24s abbrev("`v'",24) as result %8.0f `order'[1,`i'] } di as text "{hline 40}" matrix colnames `order' = `varlist' return matrix order = `order' * optionally generate the differenced (stationary) series if ("`generate'"!="") { local ords "" forvalues i = 1/`N' { local ords "`ords' `=`order'[1,`i']'" } _bu_diff `varlist' if `touse', orders(`ords') generate(`generate') replace return local newvars "`r(newvars)'" } end program define _bu_order_test, rclass * run one unit-root test and return r(rej) = 1xk rejection row-vector version 14.0 gettoken method 0 : 0 syntax varlist [if] [in], level(real) [*] local k : word count `varlist' tempname rej if ("`method'"=="ur" | ("`method'"=="adf" & `k'>1)) { if ("`method'"=="adf") { matrix `rej' = J(1,`k',.) local j = 0 foreach v of local varlist { local ++j bootur adf `v', `options' matrix `rej'[1,`j'] = (r(p_value) < `level') } } else { bootur ur `varlist', level(`level') `options' matrix P = r(p_value) matrix `rej' = J(1,`k',.) forvalues j = 1/`k' { matrix `rej'[1,`j'] = (P[`j',1] < `level') } } } else if ("`method'"=="fdr") { bootur fdr `varlist', level(`level') `options' matrix R = r(rejections) matrix `rej' = R' } else if ("`method'"=="sqt") { bootur sqt `varlist', level(`level') `options' matrix R = r(rejections) matrix `rej' = R' } else if ("`method'"=="adf") { bootur adf `varlist', `options' matrix `rej' = (r(p_value) < `level') } else if ("`method'"=="union" | "`method'"=="bootadf") { bootur `method' `varlist', level(`level') `options' matrix `rej' = (r(p_value) < `level') } else { di as error "invalid method(): `method'" exit 198 } return matrix rej = `rej' end *------------------------------------------------------------------------------- * plotmiss : visualise the pattern of missing values (plot_missing_values) *------------------------------------------------------------------------------- program define _bu_plotmiss, rclass version 14.0 syntax varlist(numeric) [if] [in], [ Name(string) TITLE(string) ] marksample touse, novarlist markout `touse', strok preserve qui keep if `touse' local N : word count `varlist' tempvar obsid qui gen `obsid' = _n local T = _N * build the long tile dataset via mata (type codes 1..4) tempname TYPE mata: bu_missmap("`varlist'", "`TYPE'") * TYPE is (T*N) x 3 : varindex, obs, typecode drop _all qui svmat double `TYPE', name(mm) rename mm1 varidx rename mm2 obs rename mm3 mtype * label the x axis with variable names local vi = 0 foreach v of local varlist { local ++vi label define _bumv `vi' "`v'", modify } label values varidx _bumv if ("`name'"=="") local name "bootur_missing" if (`"`title'"'=="") local title "Missing-value pattern" twoway /// (scatter obs varidx if mtype==1, ms(S) msize(*1.6) mcolor("77 175 74")) /// (scatter obs varidx if mtype==2, ms(S) msize(*1.6) mcolor("55 126 184")) /// (scatter obs varidx if mtype==3, ms(S) msize(*1.6) mcolor("152 78 163")) /// (scatter obs varidx if mtype==4, ms(S) msize(*1.6) mcolor("228 26 28")), /// legend(order(1 "Observed" 2 "Balanced NA" 3 "Unbalanced NA" 4 "Internal NA") /// rows(1) size(small)) /// xlabel(1(1)`N', valuelabel angle(45) noticks) /// ytitle("Observation") xtitle("") title(`"`title'"') /// graphregion(color(white)) plotregion(margin(zero)) /// name(`name', replace) restore di as text "missing-value map drawn (graph `name')" end *------------------------------------------------------------------------------- * plotorder : bar chart of the orders of integration (plot_order_integration) *------------------------------------------------------------------------------- program define _bu_plotorder, rclass version 14.0 syntax [anything] [, Name(string) TITLE(string) ] * anything: a matrix name (row vector of orders) or run after bootur order if ("`anything'"=="") local anything "r(order)" tempname M capture matrix `M' = `anything' if (_rc) { di as error "supply a (row) matrix of orders, e.g. plotorder ordmat" exit 198 } local N = colsof(`M') preserve drop _all qui set obs `N' qui gen order = . qui gen sid = _n local cn : colnames `M' forvalues i = 1/`N' { qui replace order = `M'[1,`i'] in `i' local nm : word `i' of `cn' if ("`nm'"=="") local nm "v`i'" label define _buord `i' "`nm'", modify } label values sid _buord if ("`name'"=="") local name "bootur_orders" if (`"`title'"'=="") local title "Order of integration" twoway (bar order sid, horizontal barwidth(0.7) /// fcolor("27 158 119") lcolor(white)), /// ylabel(1(1)`N', valuelabel angle(0) noticks) /// xlabel(0(1)`=max(2,`N')') xtitle("Order of integration") ytitle("") /// title(`"`title'"') graphregion(color(white)) name(`name', replace) restore di as text "orders-of-integration plot drawn (graph `name')" end version 14.0 mata: mata clear //------------------------------------------------------------------------------- // struct holding one ADF fit //------------------------------------------------------------------------------- struct bu_adf { real scalar t // ADF t-statistic real scalar c // ADF normalised (coefficient) statistic real colvector par // [gamma; lag coefficients] real colvector res // full-model residuals e (length n) real colvector bres // restricted residuals e_b (length n) } //------------------------------------------------------------------------------- // lag matrix: n x (p) with column j = x shifted down by j (first j rows 0), // then (if trim) drop the first p rows. Mirrors R lag_matrix() for a vector. //------------------------------------------------------------------------------- real matrix bu_lagmat(real colvector x, real scalar p, real scalar trim) { real scalar n, j real matrix lx n = rows(x) if (p==0) return(J(n,0,0)) lx = J(n, p, 0) for (j=0; j<=p-1; j++) { lx[|(j+2), (j+1) \ n, (j+1)|] = x[|1 \ (n-j-1)|] } if (trim) return(lx[|(p+1), 1 \ n, p|]) return(lx) } //------------------------------------------------------------------------------- // OLS coefficient vector //------------------------------------------------------------------------------- real colvector bu_ols(real colvector y, real matrix X) { return(invsym(cross(X,X)) * cross(X,y)) } //------------------------------------------------------------------------------- // diff a vector: x - c*lag(x) ; if trim drop the first observation //------------------------------------------------------------------------------- real colvector bu_diffv(real colvector x, real scalar trim, real scalar c) { real scalar n real colvector lx, dx n = rows(x) lx = J(n,1,0) lx[|2 \ n|] = x[|1 \ n-1|] dx = x - c*lx if (trim) return(dx[|2 \ n|]) return(dx) } //------------------------------------------------------------------------------- // diff a matrix (columnwise) //------------------------------------------------------------------------------- real matrix bu_diffm(real matrix x, real scalar trim, real scalar c) { real scalar n, k real matrix lx, dx n = rows(x); k = cols(x) lx = J(n,k,0) lx[|2,1 \ n,k|] = x[|1,1 \ n-1,k|] dx = x - c*lx if (trim) return(dx[|2,1 \ n,k|]) return(dx) } //------------------------------------------------------------------------------- // (two-step) detrending: OLS (QD=0) or Quasi-Difference/GLS (QD=1) // dc: 0 none, 1 intercept, 2 intercept+trend. cbar = 7 (dc1), 13.5 (dc2). //------------------------------------------------------------------------------- real colvector bu_detrend(real colvector y, real scalar dc, real scalar QD) { real scalar n, i, cn, cbar real matrix d, dQD real colvector yQD, b n = rows(y) if (dc<=0) return(y) d = J(n, dc, 0) for (i=0; i<=dc-1; i++) d[.,i+1] = (1::n):^i if (dc==1) cbar = 7 else cbar = 13.5 if (QD) cn = 1 - cbar/n else cn = 0 dQD = bu_diffm(d, 0, cn) yQD = bu_diffv(y, 0, cn) b = bu_ols(yQD, dQD) return(y - d*b) } //------------------------------------------------------------------------------- // core ADF regression (two-step detrending). Faithful port of adf_cpp(). //------------------------------------------------------------------------------- struct bu_adf scalar bu_adf_fit(real colvector z, real scalar p, real scalar dc, real scalar QD, real scalar trim, real scalar trim_ic) { real scalar n, trim_t, s2 real colvector y, ylag, ydif, b, e, eb, els real matrix X, xxi, xls real colvector yls struct bu_adf scalar out n = rows(z) y = bu_detrend(z, dc, QD) ylag = J(n,1,0); ylag[|2 \ n|] = y[|1 \ n-1|] ydif = y - ylag if (p>0) X = ylag , bu_lagmat(ydif, p, 0) else X = ylag if (trim_ic==0) trim_t = p else trim_t = trim_ic if (trim) { xls = X[|(trim_t+2),1 \ n, cols(X)|] yls = ydif[|(trim_t+2) \ n|] } else { xls = X yls = ydif } xxi = invsym(cross(xls,xls)) b = xxi * cross(xls,yls) e = ydif - X*b eb = ydif - ylag*b[1] els = e[|(trim_t+2) \ n|] s2 = cross(els,els)/rows(els) out.t = b[1]/sqrt(s2*xxi[1,1]) out.c = n*b[1]/(1 - sum(b) + b[1]) out.par = b out.res = e out.bres = eb return(out) } //------------------------------------------------------------------------------- // Nadaraya-Watson variance estimate used in the rescaled information criteria. //------------------------------------------------------------------------------- real colvector bu_npve(real colvector z, real scalar h) { real scalar n real colvector r real matrix xm, k n = rows(z) r = (1::n)/n xm = (r*J(1,n,1) - J(n,1,1)*r') / h k = normalden(xm) return( (k*(z:^2)) :/ rowsum(k) ) } real colvector bu_rescale(real colvector y, real scalar h, real scalar p, real scalar dc, real scalar QD, real scalar trim, real scalar trim_ic) { struct bu_adf scalar fit real colvector u, ydif, shat fit = bu_adf_fit(y, p, dc, QD, trim, trim_ic) u = fit.res ydif = bu_diffv(y, 0, 1) shat = sqrt(bu_npve(u, h)) return( runningsum(ydif:/shat) ) } //------------------------------------------------------------------------------- // information criterion value. ic: 1 AIC, 2 BIC, 3 MAIC, 4 MBIC. //------------------------------------------------------------------------------- real scalar bu_ic(real scalar ic, real colvector e, real scalar k, real scalar n, real scalar b0, real colvector ylag) { real scalar s2, tk, val s2 = cross(e,e)/n if (ic==1) val = log(s2) + k*2/n else if (ic==2) val = log(s2) + k*log(n)/n else { tk = cross(ylag,ylag)*(b0*b0)/s2 if (ic==3) val = log(s2) + (k+tk)*2/n else val = log(s2) + (k+tk)*log(n)/n } return(val) } // first & last non-missing row of a column (assumes no internal gaps) real rowvector bu_range(real colvector x) { real colvector idx idx = selectindex(x:!=.) return( (idx[1], idx[rows(idx)]) ) } // 1 if any series has an internal missing (a gap inside its observed span) real scalar bu_check_missing(string scalar vlist, string scalar touse) { real matrix Y real scalar N, iN, bad real colvector idx Y = st_data(., vlist, touse) N = cols(Y) bad = 0 for (iN=1; iN<=N; iN++) { idx = selectindex(Y[.,iN]:!=.) if (rows(idx) >= 2) { if ( max(idx[|2\rows(idx)|] - idx[|1\rows(idx)-1|]) != 1 ) bad = 1 } } return(bad) } // 1 if all series share the same non-missing first/last observation real scalar bu_all_balanced(string scalar vlist, string scalar touse) { real matrix Y, rng real scalar N, iN Y = st_data(., vlist, touse) N = cols(Y) rng = J(2,N,0) for (iN=1; iN<=N; iN++) rng[.,iN] = bu_range(Y[.,iN])' return( (max(rng[1,.])==min(rng[1,.])) & (max(rng[2,.])==min(rng[2,.])) ) } real scalar bu_argmin(real colvector v) { real scalar i, mi, mv mv = v[1]; mi = 1 for (i=2; i<=rows(v); i++) { if (v[i] < mv) { mv = v[i] mi = i } } return(mi) } //------------------------------------------------------------------------------- // two-step lag selection. Faithful port of adf_selectlags_cpp(). //------------------------------------------------------------------------------- struct bu_adf scalar bu_selectlags(real colvector y, real scalar pmin, real scalar pmax, real scalar ic, real scalar dc, real scalar QD, real scalar ic_scale, real scalar h_rs, real scalar trim) { real scalar n, ip, popt real colvector ys, ylag, icv struct bu_adf scalar adfp n = rows(y) ys = y if (ic_scale) ys = bu_rescale(y, h_rs, 0, dc, 0, 1, 0) ylag = bu_detrend(ys, dc, 0)[|(pmax+2) \ (n-1)|] icv = J(pmax-pmin+1, 1, 0) for (ip=pmin; ip<=pmax; ip++) { adfp = bu_adf_fit(ys, ip, dc, 0, trim, pmax) icv[ip-pmin+1] = bu_ic(ic, adfp.res[|(pmax+2) \ n|], ip, n-pmax-1, adfp.par[1], ylag) } popt = bu_argmin(icv) - 1 + pmin return( bu_adf_fit(y, popt, dc, QD, trim, 0) ) } //=============================================================================== // One-step ADF (used only by adf, two_step=FALSE) //=============================================================================== struct bu_adf scalar bu_adf1_fit(real colvector z, real scalar p, real scalar dc, real scalar trim, real scalar trim_ic) { real scalar n, i, trim_t, s2 real colvector y, ylag, ydif, b, e, els, yls real matrix X, xxi, xls, d struct bu_adf scalar out n = rows(z) y = bu_detrend(z, 0, 0) // no detrending of y ylag = J(n,1,0); ylag[|2 \ n|] = y[|1 \ n-1|] ydif = y - ylag if (p>0) X = ylag , bu_lagmat(ydif, p, 0) else X = ylag if (dc>0) { d = J(n,dc,0) for (i=0; i<=dc-1; i++) d[.,i+1] = (1::n):^i X = X , d } if (trim_ic==0) trim_t = p else trim_t = trim_ic if (trim) { xls = X[|(trim_t+2),1 \ n, cols(X)|] yls = ydif[|(trim_t+2) \ n|] } else { xls = X; yls = ydif } xxi = invsym(cross(xls,xls)) b = xxi * cross(xls,yls) e = ydif - X*b els = e[|(trim_t+2) \ n|] s2 = cross(els,els)/rows(els) out.t = b[1]/sqrt(s2*xxi[1,1]) out.c = n*b[1]/(1 - sum(b) + b[1]) out.par = b out.res = e out.bres = e return(out) } real colvector bu_rescale1(real colvector y, real scalar h, real scalar p, real scalar dc, real scalar trim, real scalar trim_ic) { struct bu_adf scalar fit real colvector u, ydif, shat fit = bu_adf1_fit(y, p, dc, trim, trim_ic) u = fit.res ydif = bu_diffv(y, 0, 1) shat = sqrt(bu_npve(u, h)) return( runningsum(ydif:/shat) ) } struct bu_adf scalar bu_selectlags1(real colvector y, real scalar pmin, real scalar pmax, real scalar ic, real scalar dc, real scalar ic_scale, real scalar h_rs, real scalar trim) { real scalar n, ip, popt real colvector ys, ylag, icv struct bu_adf scalar adfp n = rows(y) ys = y if (ic_scale) ys = bu_rescale1(y, h_rs, 0, dc, 1, 0) ylag = bu_detrend(ys, 0, 0)[|(pmax+2) \ (n-1)|] icv = J(pmax-pmin+1, 1, 0) for (ip=pmin; ip<=pmax; ip++) { adfp = bu_adf1_fit(ys, ip, dc, trim, pmax) icv[ip-pmin+1] = bu_ic(ic, adfp.res[|(pmax+2) \ n|], ip, n-pmax-1, adfp.par[1], ylag) } popt = bu_argmin(icv) - 1 + pmin return( bu_adf1_fit(y, popt, dc, trim, 0) ) } //------------------------------------------------------------------------------- // validation harness: one series, one spec //------------------------------------------------------------------------------- void _bu_check_one(string scalar var, string scalar touse, real scalar dc, real scalar ic, real scalar scale, real scalar two, real scalar pmin, real scalar pmax, real scalar qd, string scalar gnm, string scalar tnm, string scalar lnm) { real colvector y, ys, ylag, icv real scalar n, ip, pols real rowvector rg struct bu_adf scalar f, adfp y = st_data(., var, touse) rg = bu_range(y) y = y[|rg[1] \ rg[2]|] if (qd) { // union assembly: select lag under OLS, refit under QD at that lag f = bu_selectlags(y, pmin, pmax, ic, dc, 0, scale, 0.1, 1) pols = rows(f.par)-1 f = bu_adf_fit(y, pols, dc, 1, 1, 0) } else if (two) f = bu_selectlags(y, pmin, pmax, ic, dc, 0, scale, 0.1, 1) else f = bu_selectlags1(y, pmin, pmax, ic, dc, scale, 0.1, 1) st_numscalar(gnm, f.par[1]) st_numscalar(tnm, f.t) st_numscalar(lnm, rows(f.par)-1) // diagnostic: dump two-step IC vector + rescale internals if (two) { real colvector u0, ydif0, shat0 struct bu_adf scalar f0 n = rows(y) ys = y if (scale) { f0 = bu_adf_fit(y, 0, dc, 0, 1, 0) u0 = f0.res ydif0 = bu_diffv(y, 0, 1) shat0 = sqrt(bu_npve(u0, 0.1)) st_local("shmin", strofreal(colmin(shat0), "%12.6g")) st_local("umiss", strofreal(missing(u0))) ys = runningsum(ydif0:/shat0) st_local("ysmiss", strofreal(missing(ys))) } ylag = bu_detrend(ys, dc, 0)[|(pmax+2) \ (n-1)|] icv = J(pmax-pmin+1, 1, 0) for (ip=pmin; ip<=pmax; ip++) { adfp = bu_adf_fit(ys, ip, dc, 0, 1, pmax) icv[ip-pmin+1] = bu_ic(ic, adfp.res[|(pmax+2) \ n|], ip, n-pmax-1, adfp.par[1], ylag) } st_local("icvec", invtokens(strofreal(icv', "%10.5g"))) } } //=============================================================================== // BOOTSTRAP ENGINE //=============================================================================== // Observed ADF test statistics / parameter estimates / lags for all N series // over the D = |dc|*|detr| deterministic-detrending combinations, in the R // ordering: [OLS/dc1, OLS/dc2, ..., QD/dc1, QD/dc2, ...]. trim = 1. // returns nothing; fills tests (D x N), par (D x N), lags (D x N) by ref void bu_tests_panel(real matrix Y, real matrix rng, real scalar pmin, real scalar pmax, real scalar ic, real rowvector dcv, real rowvector detrv, real scalar ic_scale, real scalar h_rs, real matrix tests, real matrix par, real matrix lags) { real scalar N, nd, iN, idc, hasOLS, hasQD, p_ols, rr real colvector y struct bu_adf scalar fO, fQ real matrix tO, tQ, pO, pQ, lO N = cols(Y) nd = cols(dcv) hasOLS = anyof(detrv, 1) hasQD = anyof(detrv, 2) tO = J(nd, N, 0); pO = J(nd, N, 0); lO = J(nd, N, 0) tQ = J(nd, N, 0); pQ = J(nd, N, 0) for (iN=1; iN<=N; iN++) { y = Y[|rng[1,iN], iN \ rng[2,iN], iN|] for (idc=1; idc<=nd; idc++) { if (hasOLS) { fO = bu_selectlags(y, pmin, pmax, ic, dcv[idc], 0, ic_scale, h_rs, 1) tO[idc,iN] = fO.t pO[idc,iN] = fO.par[1] lO[idc,iN] = rows(fO.par)-1 if (hasQD) { p_ols = rows(fO.par)-1 fQ = bu_adf_fit(y, p_ols, dcv[idc], 1, 1, 0) tQ[idc,iN] = fQ.t pQ[idc,iN] = fQ.par[1] } } else { // QD only fQ = bu_selectlags(y, pmin, pmax, ic, dcv[idc], 1, ic_scale, h_rs, 1) tQ[idc,iN] = fQ.t pQ[idc,iN] = fQ.par[1] lO[idc,iN] = rows(fQ.par)-1 } } } if (hasOLS & hasQD) { tests = tO \ tQ ; par = pO \ pQ ; lags = lO \ lO } else if (hasOLS) { tests = tO ; par = pO ; lags = lO } else { tests = tQ ; par = pQ ; lags = lO } } // Observed statistics for a single bootstrap pseudo-panel y_star (D x N -> row) real rowvector bu_tests_star(real matrix Y, real matrix rng, real scalar pmin, real scalar pmax, real scalar ic, real rowvector dcv, real rowvector detrv, real scalar ic_scale, real scalar h_rs) { real matrix tests, par, lags bu_tests_panel(Y, rng, pmin, pmax, ic, dcv, detrv, ic_scale, h_rs, tests, par, lags) return( vec(tests)' ) // column-major: series-blocks of D } // DGP fit: residuals eb (restricted) and e (full), AR coefficients, per series. // Uses QD=0, trim=0, dc = dc_boot. Fills u (T x N eb), e (T x N), arest // (pmax x N) by reference; zeros outside the non-missing range. void bu_dgp_panel(real matrix Y, real matrix rng, real scalar pmin, real scalar pmax, real scalar ic, real scalar dc_boot, real scalar ic_scale, real scalar h_rs, real matrix u, real matrix e, real matrix arest) { real scalar T, N, iN, p, a, b real colvector y struct bu_adf scalar f T = rows(Y); N = cols(Y) u = J(T, N, 0); e = J(T, N, 0); arest = J(pmax, N, 0) for (iN=1; iN<=N; iN++) { a = rng[1,iN] ; b = rng[2,iN] y = Y[|a, iN \ b, iN|] f = bu_selectlags(y, pmin, pmax, ic, dc_boot, 0, ic_scale, h_rs, 0) p = rows(f.par)-1 u[|a, iN \ b, iN|] = f.bres e[|a, iN \ b, iN|] = f.res if (p>=1) arest[|1, iN \ p, iN|] = f.par[|2 \ p+1|] } } // AR(1)-driven wild sequence for AWB (returns rows(x)+1 values incl. the init) real colvector bu_genAR1(real colvector x, real scalar ar, real scalar init) { real scalar m, t real colvector y m = rows(x) y = J(m+1,1,0) y[1] = init for (t=2; t<=m+1; t++) y[t] = x[t-1] + ar*y[t-1] return(y) } // AR(p) recursion for the sieve (init zeros, x length T -> output length T) real colvector bu_genARp(real colvector x, real colvector ar) { real scalar T, p, t, j, s real colvector y, arr T = rows(x) p = rows(ar) arr = ar[p::1] // reversed y = J(T+p, 1, 0) for (t=p+1; t<=T+p; t++) { s = x[t-p] for (j=1; j<=p; j++) s = s + arr[j]*y[t-p+j-1] y[t] = s } return( y[|p+1 \ T+p|] ) } // trapezoid window and self-convolution for the DWB weight matrix real scalar bu_wtrap(real scalar x, real scalar c) { real scalar y y = 0 if (x>=0 & x=c & x<=1-c) y = 1 else if (x>1-c & x<=1) y = (1-x)/c return(y) } real scalar bu_selfconv(real scalar t, real scalar c, real scalar ngrid) { real scalar a, bnd, h, i, x, s, at at = abs(t) a = -1 ; bnd = 1 h = (bnd-a)/ngrid s = 0 for (i=0; i<=ngrid; i++) { x = a + i*h real scalar w w = bu_wtrap(x,c)*bu_wtrap(x+at,c) if (i==0 | i==ngrid) s = s + 0.5*w else s = s + w } return(s*h) } real matrix bu_dwb_s(real scalar T, real scalar l) { real scalar i, j, denom, c, ng real colvector cv, eval real matrix m, evec, s c = 0.43 ; ng = 1000 denom = bu_selfconv(0, c, ng) cv = J(2*T-1, 1, 0) // self-conv for lags -(T-1)..(T-1) for (i=1; i<=2*T-1; i++) cv[i] = bu_selfconv((i-T)/l, c, ng)/denom m = J(T,T,0) for (i=1; i<=T; i++) { for (j=1; j<=T; j++) m[i,j] = cv[(i-j)+T] } symeigensystem(m, evec, eval) eval = eval' for (i=1; i<=T; i++) { if (eval[i] <= 1e-10) eval[i] = 1e-10 } // sign-normalise eigenvectors by the sign of their first element for (j=1; j<=T; j++) { if (evec[1,j] < 0) evec[.,j] = -evec[.,j] } s = evec * diag(sqrt(eval)) return(s) } // column-wise cumulative sum of a matrix real matrix bu_colcumsum(real matrix X) { real scalar n, j real matrix Y n = rows(X) Y = X for (j=2; j<=n; j++) Y[j,.] = Y[j-1,.] + X[j,.] return(Y) } // One bootstrap pseudo-panel y_star (T x N) for a given method boot=1..6 real matrix bu_boot_dgp(real scalar boot, real matrix u, real matrix e, real colvector z, real colvector iidx, real scalar l, real matrix sD, real scalar ar, real matrix arest) { real scalar T, N, nb, iN, bnum, t real matrix ustar, ystar, xirep real colvector xi, zi, ecol T = rows(u) ; N = cols(u) if (boot==1) { // MBB nb = ceil(T/l) ustar = J(nb*l+1, N, 0) for (bnum=1; bnum<=nb; bnum++) { ustar[|(bnum-1)*l+2, 1 \ (bnum-1)*l+1+l, N|] = u[|iidx[bnum], 1 \ iidx[bnum]+l-1, N|] } ystar = bu_colcumsum(ustar) return( ystar[|rows(ystar)-T+1, 1 \ rows(ystar), N|] ) } else if (boot==2) { // BWB nb = ceil(T/l) xi = J(nb*l, 1, 0) for (bnum=1; bnum<=nb; bnum++) xi[|(bnum-1)*l+1 \ bnum*l|] = J(l,1,z[bnum]) xirep = xi[|1 \ T|] * J(1,N,1) ustar = J(1,N,0) \ (u :* xirep) ystar = bu_colcumsum(ustar) return( ystar[|2, 1 \ T+1, N|] ) } else if (boot==3) { // DWB xi = sD * z[|1 \ T|] xirep = xi * J(1,N,1) ustar = J(1,N,0) \ (u :* xirep) ystar = bu_colcumsum(ustar) return( ystar[|2, 1 \ T+1, N|] ) } else if (boot==4) { // AWB zi = z[|2 \ T|] * sqrt(1-ar*ar) xi = bu_genAR1(zi, ar, z[1]) xirep = xi * J(1,N,1) ustar = J(1,N,0) \ (u :* xirep) ystar = bu_colcumsum(ustar) return( ystar[|2, 1 \ T+1, N|] ) } else if (boot==5) { // SB real matrix estar estar = e[iidx[|1 \ T|], .] ustar = J(T, N, 0) for (iN=1; iN<=N; iN++) ustar[.,iN] = bu_genARp(estar[.,iN], arest[.,iN]) ustar = J(1,N,0) \ ustar ystar = bu_colcumsum(ustar) return( ystar[|2, 1 \ T+1, N|] ) } else { // SWB (boot==6) real matrix es2 es2 = (z * J(1,N,1)) :* e ustar = J(T, N, 0) for (iN=1; iN<=N; iN++) ustar[.,iN] = bu_genARp(es2[.,iN], arest[.,iN]) ustar = J(1,N,0) \ ustar ystar = bu_colcumsum(ustar) return( ystar[|2, 1 \ T+1, N|] ) } } // Full bootstrap loop -> t_star (B x (D*N)). rng0 is the range used for the // star tests (1..T for every column, since y_star is balanced when joint). real matrix bu_bootstrap(real scalar B, real scalar boot, real matrix u, real matrix e, real scalar l, real matrix sD, real scalar ar, real matrix arest, real scalar pmin, real scalar pmax, real scalar ic, real rowvector dcv, real rowvector detrv, real scalar ic_scale, real scalar h_rs, real matrix rng, real scalar joint, real scalar showdots) { real scalar T, N, D, ub, iB, iN, step real matrix z, iidx, tstar, ystar, rng1 real rowvector rr T = rows(u) ; N = cols(u) D = cols(dcv)*cols(detrv) if (boot==6) ub = 1 else ub = l z = rnormal(T, B, 0, 1) iidx = ceil(runiform(T, B) :* (T-ub+1)) tstar = J(B, D*N, 0) step = max((1, floor(B/40))) if (joint) { for (iB=1; iB<=B; iB++) { ystar = bu_boot_dgp(boot, u, e, z[.,iB], iidx[.,iB], l, sD, ar, arest) tstar[iB,.] = bu_tests_star(ystar, rng, pmin, pmax, ic, dcv, detrv, ic_scale, h_rs) if (showdots & mod(iB,step)==0) printf(".") } } else { for (iB=1; iB<=B; iB++) { for (iN=1; iN<=N; iN++) { ystar = bu_boot_dgp(boot, u[.,iN], e[.,iN], z[.,iB], iidx[.,iB], l, sD, ar, arest[.,iN]) rng1 = (1 \ rows(ystar)) rr = bu_tests_star(ystar, rng1, pmin, pmax, ic, dcv, detrv, ic_scale, h_rs) tstar[|iB, (iN-1)*D+1 \ iB, iN*D|] = rr } if (showdots & mod(iB,step)==0) printf(".") } } if (showdots) displayflush() return(tstar) } // individual bootstrap p-values: P*(t* < t_obs) per column real colvector bu_iadf(real rowvector testi, real matrix tstar) { real scalar M, j, B real colvector pv M = cols(testi) ; B = rows(tstar) pv = J(M,1,0) for (j=1; j<=M; j++) pv[j] = sum(tstar[.,j] :< testi[j]) / B return(pv) } // union scaling factors: prob-quantile of each series' D bootstrap columns real matrix bu_scaling(real matrix tstar, real scalar D, real scalar prob) { real scalar N, iN, B, k real matrix sc, sub real colvector col N = cols(tstar)/D ; B = rows(tstar) sc = J(D, N, 0) k = ceil(prob*B) if (k<1) k = 1 for (iN=1; iN<=N; iN++) { sub = tstar[|1,(iN-1)*D+1 \ B,iN*D|] for (col=1; col<=D; col++) { real colvector v v = sort(sub[.,col],1) sc[col,iN] = v[k] } } return(sc) } // union statistics: min over D of (-t / scaling), for a B x (D*N) matrix real matrix bu_union(real matrix t, real scalar D, real matrix s) { real scalar B, N, iB, iN real matrix ut, block B = rows(t) ; N = cols(t)/D ut = J(B, N, 0) for (iB=1; iB<=B; iB++) { for (iN=1; iN<=N; iN++) { block = t[|iB,(iN-1)*D+1 \ iB,iN*D|]' // D x 1 ut[iB,iN] = colmin( (-block) :/ s[.,iN] ) } } return(ut) } // union statistic for a single observed row (1 x D) with scaling (D x 1) real scalar bu_union1(real rowvector t, real colvector s) { return( colmin( (-t') :/ s ) ) } // k-th smallest value in a row vector real scalar bu_ksmall(real rowvector v, real scalar k) { real rowvector sv sv = sort(v', 1)' return(sv[k]) } // one BSQT step -> (p0, p1, index_j, test_j, pval_j) real rowvector bu_bsqt_step(real scalar p0, real scalar p1, real rowvector testi, real colvector ranks, real matrix tstar, real scalar N) { real scalar index_j, test_j, i, B, k, cnt real colvector ind real matrix tsub ind = sort(ranks[|p0+1 \ N|], 1) tsub = tstar[., ind] index_j = ranks[p1] test_j = testi[index_j] B = rows(tsub) k = p1 - p0 cnt = 0 for (i=1; i<=B; i++) { if (bu_ksmall(tsub[i,.], k) < test_j) cnt = cnt + 1 } return( (p0, p1, index_j, test_j, cnt/B) ) } // BSQT sequential test. Posts rej (N x 1) and the step matrix. void bu_bsqt(real colvector pvec, real rowvector testi, real matrix tstar, real scalar level, string scalar rejnm, string scalar seqnm) { real scalar N, K, phat, i, nstep real colvector ranks, rej real matrix steps N = cols(testi) ranks = order(testi', 1) K = rows(pvec)-1 phat = N nstep = K steps = J(K,5,0) for (i=1; i<=K; i++) { steps[i,.] = bu_bsqt_step(pvec[i], pvec[i+1], testi, ranks, tstar, N) if (steps[i,5] > level) { phat = steps[i,1] nstep = i i = K + 1 } } steps = steps[|1,1 \ nstep,5|] rej = J(N,1,0) if (phat > 0) rej[ranks[|1 \ phat|]] = J(phat,1,1) st_matrix(rejnm, rej) st_matrix(seqnm, steps) } // FDR step-down (Moon-Perron / Romano-Shaikh-Wolf). Faithful port of FDR_cpp. void bu_fdr(real rowvector testi, real matrix tstar, real scalar level, string scalar rejnm, string scalar seqnm) { real scalar N, B, j, no_r, phat, pos, r, a, cc, pseudo, mrow real colvector ranks, cvfdr, rej, FDRest, cum, cmp, rj, rc1 real matrix tsub, sorted, cvstar, norej, norej_jp, seq N = cols(testi) B = rows(tstar) ranks = order(testi', 1) pseudo = 1e300 cvfdr = J(N,1,0) for (j=0; j<=N-1; j++) { tsub = tstar[., ranks[|N-j \ N|]] // B x (j+1) sorted = J(B, j+1, 0) for (r=1; r<=B; r++) sorted[r,.] = sort(tsub[r,.]', 1)' rc1 = order(sorted[.,1], 1) cvstar = sorted[rc1,.] \ J(1, j+1, pseudo) // (B+1) x (j+1) if (j==0) { pos = min((B, floor(N*level*B))) if (N*level <= 1) pos = pos - 1 cvfdr[N] = cvstar[pos+1, 1] // C++ 0-based -> +1 } else if (j < N-1) { norej_jp = J(B+1, j, 1) for (a=1; a<=B+1; a++) { for (cc=1; cc<=j; cc++) { if (cvstar[a, cc+1] <= cvfdr[N-j+cc]) norej_jp[a,cc] = 0 } } norej_jp = norej_jp , J(B+1,1,1) norej = J(B+1, j+1, 1) norej[.,1] = norej_jp[.,1] FDRest = bu_rowprod(norej) / (N - j) for (no_r=1; no_r<=j; no_r++) { norej[.,no_r] = J(B+1,1,1) - norej[.,no_r] norej[.,no_r+1] = norej_jp[.,no_r+1] FDRest = FDRest + (no_r+1) * bu_rowprod(norej) / (N - j + no_r) } cum = runningsum(FDRest / B) pos = sum(cum :<= level) cvfdr[N-j] = cvstar[pos, 1] } else { cvfdr[1] = cvstar[floor(level*B), 1] } } // rejections: cumulative product of (test_i(ranks) < cvfdr), scattered back cmp = (testi[ranks]' :< cvfdr) rj = bu_cumprod(cmp) rej = J(N,1,0) rej[ranks] = rj phat = sum(rej) mrow = min((N, phat+1)) seq = J(mrow, 3, 0) for (j=1; j<=mrow; j++) { seq[j,1] = ranks[j] seq[j,2] = testi[ranks[j]] seq[j,3] = cvfdr[j] } st_matrix(rejnm, rej) st_matrix(seqnm, seq) } real colvector bu_rowprod(real matrix X) { real scalar n, m, i, j real colvector p n = rows(X) ; m = cols(X) p = J(n,1,1) for (i=1; i<=n; i++) { for (j=1; j<=m; j++) p[i] = p[i]*X[i,j] } return(p) } real colvector bu_cumprod(real colvector v) { real scalar n, i real colvector p n = rows(v) p = J(n,1,0) p[1] = v[1] for (i=2; i<=n; i++) p[i] = p[i-1]*v[i] return(p) } //=============================================================================== // MAIN DRIVER for boot_ur / boot_fdr / boot_sqt / boot_panel //=============================================================================== // opt = (boot,B,l,arAWB,ic,ic_scale,h_rs,pmin,pmax,union,dc,detr,level,joint,mode) // mode: 1=ur, 2=fdr, 3=sqt, 4=panel void bu_run(string scalar vlist, string scalar touse, real rowvector opt, string scalar pvecnm) { real scalar boot, B, l, arAWB, ic, ic_scale, h_rs, pmin, pmax, uni real scalar dc1, detr1, level, joint, mode, T, N, D, dc_boot, iN, showdots real matrix Y, rng, tests, par, lags, tstar, sD, u, e, arest, ustar real matrix sc, pind real rowvector dcv, detrv, uobs real colvector pv, est, stat, pvalc, gm_star real scalar gm, gmp boot=opt[1]; B=opt[2]; l=opt[3]; arAWB=opt[4]; ic=opt[5]; ic_scale=opt[6] h_rs=opt[7]; pmin=opt[8]; pmax=opt[9]; uni=opt[10]; dc1=opt[11]; detr1=opt[12] level=opt[13]; joint=opt[14]; mode=opt[15]; showdots=opt[16] Y = st_data(., vlist, touse) T = rows(Y); N = cols(Y) rng = J(2,N,0) for (iN=1; iN<=N; iN++) rng[.,iN] = bu_range(Y[.,iN])' if (uni) { dcv = (1,2); detrv = (1,2); dc_boot = 2 } else { dcv = (dc1); detrv = (detr1); dc_boot = dc1 } D = cols(dcv)*cols(detrv) if (boot==3) sD = bu_dwb_s(T, l) else sD = J(1,1,0) bu_tests_panel(Y, rng, pmin, pmax, ic, dcv, detrv, ic_scale, h_rs, tests, par, lags) bu_dgp_panel(Y, rng, pmin, pmax, ic, dc_boot, ic_scale, h_rs, u, e, arest) tstar = bu_bootstrap(B, boot, u, e, l, sD, arAWB, arest, pmin, pmax, ic, dcv, detrv, ic_scale, h_rs, rng, joint, showdots) // individual (per-combination) p-values, reshaped to N x D pv = bu_iadf(vec(tests)', tstar) pind = rowshape(pv, N) // observed statistic & estimate to report at series level est = J(N,1,.) stat = J(N,1,0) pvalc = J(N,1,.) if (uni) { sc = bu_scaling(tstar, D, level) ustar = bu_union(tstar, D, sc) // B x N uobs = J(1,N,0) for (iN=1; iN<=N; iN++) uobs[iN] = bu_union1(tests[.,iN]', sc[.,iN]) stat = uobs' // estimate stays missing for union } else { stat = tests[1,.]' est = par[1,.]' } if (mode==1) { // ur if (uni) pvalc = bu_iadf(uobs, ustar) else pvalc = bu_iadf(tests[1,.], tstar) } else if (mode==2) { // fdr if (uni) bu_fdr(uobs, ustar, level, "__bu_rej", "__bu_seq") else bu_fdr(tests[1,.], tstar, level, "__bu_rej", "__bu_seq") } else if (mode==3) { // sqt real colvector pvec pvec = st_matrix(pvecnm)' if (uni) bu_bsqt(pvec, uobs, ustar, level, "__bu_rej", "__bu_seq") else bu_bsqt(pvec, tests[1,.], tstar, level, "__bu_rej", "__bu_seq") } else { // panel group-mean if (uni) { gm = mean(uobs') gm_star = bu_rowmean(ustar) } else { gm = mean(tests[1,.]') gm_star = bu_rowmean(tstar) } gmp = sum(gm_star :< gm)/B st_numscalar("__bu_gm", gm) st_numscalar("__bu_gmp", gmp) } // post everything for the ado to format st_matrix("__bu_est", est) st_matrix("__bu_stat", stat) st_matrix("__bu_pval", pvalc) st_matrix("__bu_dstat", tests') // N x D individual statistics st_matrix("__bu_dpar", par') // N x D individual estimates st_matrix("__bu_dlag", lags') // N x D individual lags st_matrix("__bu_dpval", pind) // N x D individual p-values st_numscalar("__bu_D", D) st_numscalar("__bu_N", N) st_numscalar("__bu_jointused", joint) } real colvector bu_rowmean(real matrix X) { return( rowsum(X)/cols(X) ) } // difference a column d times (row order), leaving d leading missings real colvector bu_ndiff(real colvector y, real scalar d) { real scalar k, t, n real colvector x, nx x = y for (k=1; k<=d; k++) { n = rows(x) nx = J(n,1,.) for (t=2; t<=n; t++) nx[t] = x[t] - x[t-1] x = nx } return(x) } void bu_diff_store(string scalar vin, string scalar vout, string scalar touse, real scalar d) { real colvector y, dd, idx y = st_data(., vin, touse) dd = bu_ndiff(y, d) idx = selectindex(st_data(., touse):!=0) st_store(idx, st_varindex(vout), dd) } // build a missing-value classification map (types 1..4) as (T*N) x 3 void bu_missmap(string scalar vlist, string scalar matname) { real matrix Y, rng, out real scalar T, N, iN, t, minidx, maxidx, r, ty Y = st_data(., vlist) T = rows(Y); N = cols(Y) rng = J(2,N,0) for (iN=1; iN<=N; iN++) rng[.,iN] = bu_range(Y[.,iN])' minidx = min(rng[1,.]); maxidx = max(rng[2,.]) out = J(T*N, 3, 0) r = 0 for (iN=1; iN<=N; iN++) { for (t=1; t<=T; t++) { r = r + 1 out[r,1] = iN out[r,2] = T - t + 1 // so obs 1 plots at the top ty = 1 if (t < minidx | t > maxidx) ty = 2 else if (Y[t,iN]==.) { if (t > rng[1,iN] & t < rng[2,iN]) ty = 4 else ty = 3 } out[r,3] = ty } } st_matrix(matname, out) } // asymptotic ADF for a single series -> gamma, tstat, lag, effective TT void bu_adf_run(string scalar vlist, string scalar touse, real scalar dc, real scalar ic, real scalar scale, real scalar two, real scalar pmin, real scalar pmax, string scalar gnm, string scalar tnm, string scalar lnm, string scalar ttnm) { real colvector y real rowvector rg struct bu_adf scalar f y = st_data(., vlist, touse) rg = bu_range(y) y = y[|rg[1] \ rg[2]|] if (two) f = bu_selectlags(y, pmin, pmax, ic, dc, 0, scale, 0.1, 1) else f = bu_selectlags1(y, pmin, pmax, ic, dc, scale, 0.1, 1) st_numscalar(gnm, f.par[1]) st_numscalar(tnm, f.t) st_numscalar(lnm, rows(f.par)-1) st_numscalar(ttnm, rows(y)) } // <<>> // ==== MacKinnon (1996) unit-root p-value response surfaces (from urca) ==== // (raw Mata snippet, injected into bootur.ado's engine block) real matrix bu_urc_nc() { string scalar s real rowvector v s = "" s = s + " -3.8929681 -10.711437 -55.619955 0.0033188192 -3.7185872 -9.5896921 -34.166532 0.0024951489 -3.4759164 -7.8746051 -17.490657 0.0015559156 -3.2842708 -6.4108912 -11.983755 0.0011445964 -3.0830457 -5.0567633 -6.8847143 0.00084100765 -2.9598602 -4.2584716 -5.8142463 0.00071049934 -2.8698825 -3.728955 -5.8714971 0.00062391151 -2.7977077 -3.4146427 -3.76153 0.00058693244 -2.738005 -3.1208498 -3.0462883 0.00052576407 -2.6868436 -2.8722476 -2.5438563 0.00049052976" s = s + " -2.6416852 -2.6793402 -1.8769783 0.00046804316 -2.6012597 -2.5000004 -1.6307933 0.00044346134 -2.5649446 -2.3127399 -2.0306605 0.00043039003 -2.4206178 -1.7291886 -1.3466493 0.00036467475 -2.3133304 -1.3423615 -1.3432257 0.00032107164 -2.2273143 -1.0746264 -0.8762427 0.00029789755 -2.1549047 -0.84657628 -1.0780188 0.00027830583 -2.0921655 -0.67744316 -1.0430697 0.00026027602 -2.0364584 -0.54192487 -0.96965194 0.00025249134 -1.9864146 -0.4161689 -1.1334045 0.00024079995" s = s + " -1.9407684 -0.33025231 -0.67821426 0.00023323022 -1.8988585 -0.22970773 -0.85640849 0.00022075141 -1.8598488 -0.14639728 -1.0428806 0.00021334955 -1.8233573 -0.089663803 -0.80855945 0.0002096301 -1.7891398 -0.036466865 -0.63243896 0.00020834232 -1.7567993 0.011402654 -0.51066497 0.00020539933 -1.7262105 0.066220956 -0.63708003 0.00019964762 -1.6969838 0.10236824 -0.46672787 0.00019574132 -1.6690859 0.14265332 -0.45061418 0.00018888062 -1.6424271 0.18701823 -0.66084335 0.00018426555" s = s + " -1.6167279 0.21626789 -0.64616991 0.00018128089 -1.592052 0.24778988 -0.63597502 0.00017954902 -1.5682114 0.27323562 -0.61618157 0.00017773495 -1.545274 0.29991006 -0.60544565 0.00017142168 -1.5230643 0.33148333 -0.76580637 0.00016769202 -1.5014865 0.35136387 -0.69024691 0.00016887147 -1.48057 0.37478667 -0.72930388 0.00016872218 -1.4602522 0.39420657 -0.70443601 0.00016401731 -1.4404242 0.41028681 -0.64338518 0.00016224607 -1.4210785 0.421192 -0.555484 0.00016193196" s = s + " -1.4022406 0.43747746 -0.5767082 0.00016097849 -1.3837876 0.44575298 -0.43236279 0.00015803832 -1.3657717 0.45528883 -0.35165569 0.00015694908 -1.3481032 0.46306181 -0.31093184 0.00015603528 -1.3309144 0.47972473 -0.3912126 0.00015557033 -1.3139547 0.48814175 -0.38006723 0.00015236727 -1.2973931 0.49859437 -0.35456698 0.00014909593 -1.281103 0.50420631 -0.2558317 0.0001479347 -1.2651262 0.51429788 -0.28893557 0.00014655593 -1.2494004 0.52135123 -0.27882309 0.00014412602" s = s + " -1.2339522 0.52756595 -0.23217506 0.00014437838 -1.2187168 0.53480958 -0.26406268 0.00014432408 -1.2037702 0.54338564 -0.29729024 0.00014262955 -1.1890155 0.54817266 -0.24900712 0.00014111421 -1.1744677 0.55263933 -0.22129035 0.00014067648 -1.1601444 0.55919565 -0.23623151 0.00013973685 -1.1460177 0.56517698 -0.24565447 0.00013926384 -1.1320096 0.56603215 -0.17524226 0.0001396978 -1.1182292 0.57011095 -0.19439876 0.00013916343 -1.1045866 0.57209736 -0.16074907 0.00013818167" s = s + " -1.0911074 0.5723485 -0.10524258 0.00013722018 -1.0777556 0.57229676 -0.032550732 0.00013619288 -1.0645684 0.57302396 0.014379749 0.0001355341 -1.0515244 0.5736259 0.08463918 0.00013526527 -1.0386317 0.57725998 0.07359064 0.00013519196 -1.0258492 0.5804434 0.051060038 0.00013379316 -1.01317 0.58385933 0.0061245025 0.0001342836 -1.000627 0.58959535 -0.064975446 0.00013403253 -0.9881634 0.58962904 -0.03896917 0.00013388645 -0.97580892 0.59228543 -0.092028608 0.00013256283" s = s + " -0.96353988 0.59023118 -0.025245537 0.00013350108 -0.95135437 0.59235081 -0.066927698 0.00013166673 -0.93927127 0.59495523 -0.11772564 0.00013109874 -0.92725721 0.59465251 -0.098842151 0.00013075847 -0.91538294 0.59862948 -0.1337387 0.00012991997 -0.90351591 0.59736071 -0.11575297 0.00013004651 -0.8917174 0.59827558 -0.11376391 0.00012981779 -0.87999358 0.59854644 -0.11262976 0.0001293793 -0.86836393 0.60125923 -0.1570054 0.00012809482 -0.8567663 0.60235243 -0.15341154 0.0001272749" s = s + " -0.84521821 0.60724747 -0.28142705 0.00012854899 -0.83368446 0.60414764 -0.22759417 0.00012886527 -0.82217754 0.5988247 -0.1036938 0.0001272911 -0.81077477 0.60266754 -0.1923692 0.00012744791 -0.79938013 0.60203622 -0.1545566 0.00012851757 -0.7879978 0.60310963 -0.21236507 0.0001284728 -0.77666556 0.60512141 -0.24195949 0.00012867783 -0.76531152 0.6016104 -0.18435951 0.00012764627 -0.75400509 0.60005644 -0.15758711 0.00012719162 -0.7427255 0.60263001 -0.19073741 0.00012737733" s = s + " -0.7314095 0.60169755 -0.1680994 0.00012901165 -0.72010821 0.60073899 -0.17118517 0.00012989901 -0.70881438 0.60281341 -0.19734629 0.00013096756 -0.69752254 0.605741 -0.2390057 0.00013200027 -0.68621242 0.60709003 -0.25007855 0.00013205742 -0.67486899 0.60595186 -0.21780355 0.00013308173 -0.6635531 0.60879991 -0.23907614 0.00013334734 -0.65217906 0.6093449 -0.22402739 0.00013313113 -0.6407847 0.61045036 -0.22525697 0.00013423481 -0.62935668 0.61165455 -0.18555733 0.00013561164" s = s + " -0.61784992 0.60906224 -0.1021716 0.00013524395 -0.60634066 0.61133645 -0.093563312 0.00013574028 -0.59479908 0.61243094 -0.018009253 0.00013677974 -0.58317131 0.61383424 0.033976678 0.00013851627 -0.57147764 0.61383949 0.12970947 0.0001391928 -0.55974664 0.61815202 0.13567448 0.00013865735 -0.54796415 0.62425007 0.12681084 0.00014027742 -0.53605781 0.62788652 0.1392266 0.00014160362 -0.52407823 0.63003047 0.2001274 0.00014263868 -0.51203077 0.63456809 0.22849852 0.00014476471" s = s + " -0.49993482 0.6410294 0.23651226 0.00014745515 -0.48774646 0.64854569 0.1918391 0.00014882984 -0.47546258 0.65391229 0.20689462 0.0001506309 -0.46312562 0.6619045 0.15084854 0.00015052565 -0.45068236 0.67032467 0.10052172 0.00015261383 -0.43816571 0.67906278 0.018043158 0.00015193656 -0.42553012 0.68596166 -0.034373075 0.00015324255 -0.41278915 0.68751627 0.0018499262 0.00015211374 -0.39996901 0.69262508 -0.044106098 0.00015262241 -0.38705587 0.69470103 -0.03722354 0.00015233648" s = s + " -0.37407697 0.69915156 -0.068104682 0.00015310849 -0.361036 0.70418582 -0.098476385 0.0001541773 -0.34792319 0.70730664 -0.10085193 0.00015510072 -0.33471267 0.71104064 -0.12671226 0.00015539396 -0.32142107 0.71478513 -0.15318544 0.00015678198 -0.30803148 0.71372891 -0.10938638 0.0001588677 -0.29463058 0.72056282 -0.22811198 0.00015869168 -0.28109052 0.72106303 -0.20822006 0.00015986958 -0.26754399 0.73003235 -0.37500115 0.00016005206 -0.25384762 0.7308834 -0.36414154 0.00015982446" s = s + " -0.24009915 0.73329574 -0.39708747 0.00016064215 -0.22627489 0.735505 -0.38861738 0.00016028259 -0.21236427 0.73562685 -0.37745695 0.00015906987 -0.19835027 0.73617467 -0.37412904 0.00015881342 -0.18428951 0.73644506 -0.34303983 0.00016010471 -0.17011954 0.73618375 -0.29890353 0.00016091489 -0.1558682 0.73373246 -0.2270809 0.00016122387 -0.14152198 0.73271165 -0.21753843 0.00015985951 -0.12710086 0.73048568 -0.16997032 0.00015998493 -0.11258798 0.72754793 -0.10077334 0.00015998672" s = s + " -0.09797607 0.72884778 -0.1089758 0.00016140512 -0.083287262 0.73024638 -0.10172413 0.00016110431 -0.06849948 0.73431662 -0.1856636 0.00016282783 -0.053596308 0.73398093 -0.15607347 0.00016572782 -0.03859197 0.73389218 -0.17541939 0.00016622391 -0.023470382 0.73281912 -0.18774382 0.00016697962 -0.0082414051 0.73243552 -0.16681816 0.00016735938 0.0070938188 0.73311643 -0.16926117 0.00017001641 0.022558531 0.73323536 -0.16912639 0.0001713441 0.03812404 0.73423643 -0.21631637 0.00017166919" s = s + " 0.053822876 0.73695639 -0.24374285 0.00017135216 0.069723974 0.73497901 -0.21283984 0.0001717885 0.085748648 0.73353054 -0.22887618 0.00017378309 0.10194776 0.72751724 -0.14073615 0.00017367948 0.11822374 0.72919921 -0.19181418 0.00017409813 0.13466118 0.72841873 -0.20592264 0.00017404024 0.15123202 0.72748214 -0.142825 0.00017445653 0.16799709 0.72206466 -0.015965775 0.00017698924 0.1849277 0.71874017 0.066004976 0.00017778151 0.20202274 0.71934368 0.025624359 0.00017955598" s = s + " 0.21927304 0.72039275 0.02558038 0.00018184571 0.23678102 0.71748717 0.094060252 0.00018364653 0.25451636 0.70885897 0.25088485 0.00018473165 0.27238855 0.70503684 0.37460085 0.00018585039 0.2904553 0.70469131 0.43622979 0.00018669006 0.30876289 0.69991385 0.57078447 0.0001873182 0.32730986 0.69817689 0.65322255 0.00018881434 0.34608201 0.6995444 0.65207222 0.00018932827 0.3650916 0.7011835 0.66744563 0.00019036477 0.3843845 0.70479523 0.63848692 0.00019104889" s = s + " 0.40399054 0.70127591 0.77381319 0.00019254653 0.42381859 0.70578585 0.77661271 0.00019357949 0.44400981 0.70845188 0.76564843 0.000192708 0.46452123 0.71504692 0.66984153 0.00019399621 0.48536196 0.71727201 0.7293499 0.00019557474 0.50650098 0.73110911 0.54567534 0.00019837528 0.52811782 0.72881888 0.7389034 0.00019966728 0.55006427 0.73523577 0.74836455 0.00019991704 0.57244866 0.74134438 0.80091013 0.00019902737 0.59531037 0.74367625 0.91978665 0.00020239438" s = s + " 0.61869583 0.74658014 0.99643135 0.00020585053 0.6424973 0.76305579 0.87771089 0.00020526129 0.66691707 0.76996444 0.97362052 0.00020866944 0.69183128 0.78557766 0.88754062 0.00020955134 0.71744906 0.80040433 0.80308001 0.00021295825 0.74377276 0.8121428 0.83723176 0.00021328053 0.77080898 0.82540821 0.85895798 0.00021879241 0.79865774 0.83593059 0.91965289 0.00022425873 0.8274514 0.84399837 1.0515597 0.00022793382 0.8570601 0.8606473 1.1167888 0.00023196873" s = s + " 0.88773648 0.8764175 1.1796358 0.000239097 0.91952791 0.89954267 1.1230715 0.00023965676 0.95259765 0.91572736 1.3446001 0.00024244577 0.98685842 0.94668499 1.3083248 0.00024366251 1.0226723 0.97271007 1.525428 0.00024879761 1.0602764 1.0038143 1.5791882 0.00025394977 1.0998111 1.0352814 1.8038552 0.00025746326 1.1415718 1.0682334 2.0945465 0.00025962231 1.1857427 1.1068815 2.4142752 0.0002629678 1.2329803 1.1470414 2.8497297 0.0002668687" s = s + " 1.2836101 1.2111123 2.9446612 0.000274748 1.3383168 1.2755673 3.3808455 0.00028272951 1.3981251 1.3492646 3.800062 0.00029185651 1.4640642 1.4544896 4.1103714 0.00031385459 1.5383455 1.5713643 4.61051 0.00032810887 1.6233569 1.7314438 5.1765144 0.00035016457 1.7241804 1.926988 6.350803 0.00038437873 1.8485834 2.2458904 7.2317207 0.00041561758 2.0150248 2.7264767 8.8446031 0.00048601732 2.0567136 2.8527315 9.6184847 0.00051796899" s = s + " 2.1027347 3.0008151 10.382866 0.00053898121 2.1540126 3.1822666 11.07189 0.0005577211 2.2124108 3.345449 12.865521 0.00058784437 2.279742 3.6015037 14.01974 0.00063354539 2.3603355 3.9707581 14.523069 0.00071569175 2.4616852 4.4258179 15.562243 0.00077596696 2.6009217 4.9197115 22.640137 0.00090787028 2.8233497 6.3045355 26.417722 0.0011716504 3.0328032 7.5631954 35.087491 0.0016240091 3.2991405 9.4748927 53.382413 0.0023643422" s = s + " 3.4781086 10.960181 65.777105 0.0031382626" v = strtoreal(tokens(s)) return(rowshape(v, 221)) } real matrix bu_urc_c() { string scalar s real rowvector v s = "" s = s + " -4.6498737 -19.585128 -134.04859 0.0028987599 -4.4931648 -16.807926 -122.0471 0.0022083 -4.2675648 -14.425831 -87.644297 0.0014140247 -4.0912176 -12.441367 -68.492333 0.0010882695 -3.9058347 -10.513102 -51.918428 0.0007795883 -3.7927154 -9.4090272 -44.208442 0.00066318864 -3.7100105 -8.6780936 -37.830866 0.00059497541 -3.643575 -8.1822049 -32.355767 0.00055007457 -3.589081 -7.6729277 -30.421975 0.00050967313 -3.5421556 -7.3119578 -27.319567 0.00048484389" s = s + " -3.5007351 -6.9999449 -25.038104 0.00045285603 -3.4637847 -6.7024239 -23.833155 0.00043212199 -3.4301855 -6.4601835 -22.150622 0.0004062553 -3.2980044 -5.5041761 -17.079587 0.00035316546 -3.2003475 -4.8184935 -14.518699 0.00031144263 -3.1220215 -4.3156408 -12.363387 0.00028433742 -3.0560484 -3.9360126 -10.355988 0.00027365317 -2.998901 -3.6004495 -9.2929232 0.00025607319 -2.9484715 -3.3101354 -8.2865912 0.00024419727 -2.9028807 -3.0689608 -7.3753301 0.00023798936" s = s + " -2.8613794 -2.858389 -6.5356937 0.0002256829 -2.8231707 -2.6704867 -5.8684031 0.00022093751 -2.787681 -2.5009584 -5.2732527 0.00021380944 -2.7545889 -2.3414925 -4.8579683 0.00020908348 -2.723494 -2.1950703 -4.5079621 0.00020412734 -2.6940314 -2.0720578 -4.0170328 0.00019962789 -2.6661673 -1.9520849 -3.6714342 0.00019401044 -2.6395787 -1.8455989 -3.2400536 0.00019334773 -2.6142084 -1.7412463 -2.9756301 0.00018854058 -2.589969 -1.6382011 -2.818905 0.00018342614" s = s + " -2.5666852 -1.5407058 -2.6824738 0.00017867718 -2.5442992 -1.4455909 -2.6141252 0.00017529619 -2.5226529 -1.3636472 -2.4331178 0.00017198577 -2.5017734 -1.2860456 -2.2376568 0.00016710155 -2.4815082 -1.214656 -2.002867 0.00016209084 -2.4618127 -1.1438367 -1.9410583 0.0001581832 -2.4427398 -1.0795815 -1.7379897 0.00015707067 -2.4242373 -1.0129871 -1.6157265 0.00015540008 -2.4061235 -0.95863004 -1.3584182 0.00015496957 -2.3884705 -0.90317855 -1.1645286 0.00015380726" s = s + " -2.3712834 -0.8489709 -1.0361953 0.00015197444 -2.3544588 -0.79871409 -0.89103868 0.00015155235 -2.3380393 -0.7470885 -0.82182554 0.00015056603 -2.3219515 -0.70155967 -0.63796928 0.00014928147 -2.3061761 -0.65942293 -0.43314474 0.00014830706 -2.290754 -0.6128591 -0.36048529 0.0001472114 -2.2756297 -0.56779943 -0.29121491 0.00014719167 -2.260788 -0.52409236 -0.25606839 0.00014580297 -2.2462575 -0.47609723 -0.30994831 0.00014528649 -2.2318566 -0.43929292 -0.22440805 0.00014505919" s = s + " -2.217773 -0.39701536 -0.24856318 0.00014188639 -2.2038352 -0.36444851 -0.11354554 0.0001417164 -2.1901914 -0.32532352 -0.093216784 0.00014218242 -2.1767391 -0.28747114 -0.10310017 0.00014155611 -2.1634883 -0.25016662 -0.10675131 0.00014050226 -2.1503874 -0.21601205 -0.10476651 0.00014015033 -2.1374432 -0.1858336 -0.051211317 0.00014011318 -2.124642 -0.15683431 -0.029204874 0.00013883378 -2.1120088 -0.13213346 0.12852807 0.00013857794 -2.0995474 -0.099815891 0.099058341 0.00013742797" s = s + " -2.0872295 -0.070378462 0.10800415 0.00013793884 -2.0750609 -0.042044889 0.14101216 0.00013621927 -2.0630352 -0.015141684 0.18865304 0.00013639356 -2.051136 0.015952055 0.10012626 0.00013490015 -2.0393918 0.044522529 0.086565146 0.00013472776 -2.0277077 0.067059234 0.16819054 0.00013407648 -2.016149 0.090904806 0.19550829 0.00013409339 -2.004692 0.11665323 0.15704451 0.00013291655 -1.9933276 0.13968718 0.17247884 0.00013151197 -1.9820619 0.16173157 0.19898819 0.00013162402" s = s + " -1.970886 0.18237481 0.25236052 0.00013061907 -1.9598292 0.20391525 0.28727209 0.00013018764 -1.9488223 0.22372134 0.32812856 0.00012990573 -1.9379266 0.24774843 0.26125313 0.0001296202 -1.9270853 0.26745749 0.2697937 0.00012787068 -1.9163151 0.28673153 0.28767581 0.0001281965 -1.9056274 0.30480271 0.31440368 0.00012750112 -1.8950092 0.3248774 0.2928928 0.00012633421 -1.8844373 0.34183301 0.33094959 0.0001261714 -1.8739396 0.36093334 0.29415769 0.00012542124" s = s + " -1.8635243 0.38142067 0.22614552 0.00012490566 -1.8531302 0.39716105 0.24303066 0.0001248705 -1.8427987 0.4096043 0.36111823 0.00012450876 -1.832549 0.42827858 0.35044668 0.00012429566 -1.822312 0.44312533 0.38373089 0.00012339286 -1.8121476 0.45825395 0.41736017 0.00012354584 -1.8020189 0.47456176 0.39780568 0.00012299378 -1.7919441 0.49156607 0.37160352 0.00012115398 -1.7818927 0.50669755 0.35931258 0.00012090362 -1.7718547 0.52020078 0.3685279 0.00011990979" s = s + " -1.7618784 0.53552131 0.37399927 0.00012096096 -1.7519035 0.5466972 0.42995618 0.00012096099 -1.7419868 0.56180886 0.40902997 0.00012127621 -1.7321082 0.57467469 0.42868053 0.0001203501 -1.7222375 0.58708826 0.45175084 0.00012043651 -1.7123798 0.59917697 0.48860303 0.00012051206 -1.7025509 0.61170671 0.49474275 0.00012071179 -1.6927408 0.62506026 0.46000167 0.00012050594 -1.6829693 0.63659023 0.49564593 0.00012111163 -1.6731939 0.64805654 0.50804098 0.00012038583" s = s + " -1.6634091 0.65672036 0.55607608 0.00012102541 -1.6536497 0.66769315 0.55911498 0.00012149389 -1.6439062 0.67742799 0.59076277 0.00012261731 -1.6341886 0.68994994 0.56825647 0.00012213119 -1.6244687 0.70100834 0.57240705 0.00012293942 -1.6147364 0.71044565 0.60084676 0.00012159089 -1.6050132 0.72047181 0.61499829 0.00012113944 -1.5952755 0.72923289 0.6508903 0.00012120266 -1.5855407 0.7394689 0.6451626 0.00012082145 -1.575834 0.7491741 0.64206303 0.00012086709" s = s + " -1.5661226 0.76123129 0.60141833 0.00012178366 -1.5563725 0.76883728 0.64611202 0.00012207516 -1.5466262 0.77749522 0.67164232 0.00012318791 -1.5368887 0.79044802 0.59971309 0.0001234695 -1.5271193 0.80044531 0.60613339 0.00012451886 -1.5173571 0.80941987 0.6172132 0.00012384741 -1.5075829 0.82011412 0.59323485 0.00012382238 -1.4977672 0.83044274 0.58712998 0.0001238026 -1.4879122 0.83904973 0.59087687 0.00012307983 -1.4780499 0.84732674 0.61608062 0.0001226664" s = s + " -1.4681619 0.8567715 0.60979842 0.00012252223 -1.4582463 0.86435803 0.64988385 0.00012330911 -1.4482883 0.87451774 0.61543666 0.00012347216 -1.4382813 0.88263113 0.62973992 0.00012380557 -1.4282536 0.88956398 0.67468172 0.00012417842 -1.4182008 0.89929049 0.67428793 0.00012422941 -1.4080885 0.90829942 0.67159819 0.00012506386 -1.3979342 0.91703848 0.66877596 0.00012625226 -1.3877694 0.92859017 0.61995067 0.00012699355 -1.377523 0.93604289 0.63846649 0.00012666064" s = s + " -1.3672057 0.94426847 0.6251151 0.00012737116 -1.3568339 0.95188127 0.63766879 0.00012721756 -1.3464079 0.95959007 0.66127287 0.00012944453 -1.3359058 0.96366545 0.74578117 0.000130002 -1.3253685 0.97110447 0.79643193 0.00013007992 -1.3147624 0.97855405 0.83213055 0.00013042507 -1.3041106 0.98835692 0.81911505 0.00013139037 -1.2933783 0.99697287 0.86808312 0.00013181101 -1.2825497 1.0047337 0.89757564 0.00013204824 -1.2716606 1.0146399 0.90814494 0.00013257959" s = s + " -1.2606975 1.0258937 0.87659985 0.00013332493 -1.249625 1.0349109 0.88896456 0.00013332758 -1.2384725 1.0438386 0.91800985 0.00013398711 -1.2271853 1.0534207 0.9059538 0.0001346633 -1.2158071 1.0597995 1.0029103 0.00013481324 -1.2043342 1.0694094 1.0105287 0.00013585221 -1.1927374 1.0817023 0.96515325 0.00013634858 -1.1810471 1.0916938 1.0085506 0.00013694656 -1.1692236 1.1026291 1.026851 0.00013788806 -1.1572565 1.1113878 1.063002 0.00013883535" s = s + " -1.1451871 1.1272707 0.97342183 0.00013931866 -1.1329554 1.1390215 0.97571997 0.00014041549 -1.1205568 1.1488008 1.0272798 0.0001409265 -1.1080251 1.1621514 0.9981148 0.00014220521 -1.0953689 1.1751472 0.99909208 0.00014229817 -1.0825027 1.184276 1.1029791 0.00014424766 -1.0694471 1.1955497 1.1515205 0.00014635427 -1.0562396 1.2073469 1.1921759 0.00014761855 -1.0428562 1.223239 1.1346513 0.00014852001 -1.029262 1.2348075 1.1945083 0.00014902807" s = s + " -1.0154635 1.2433879 1.3252913 0.00014947863 -1.0014375 1.2567831 1.3560022 0.00015208154 -0.98722361 1.2729485 1.3623659 0.00015296326 -0.97276525 1.2849843 1.4498376 0.000153746 -0.95804968 1.2990845 1.4649926 0.00015448352 -0.94307789 1.3116712 1.5199257 0.00015402427 -0.92788147 1.3271247 1.5335271 0.00015443462 -0.91240062 1.3458079 1.4685748 0.00015632675 -0.89662363 1.3603465 1.5140667 0.0001568523 -0.88051439 1.3707204 1.6170701 0.00015674102" s = s + " -0.86415213 1.390553 1.5259264 0.00015760527 -0.84741271 1.4046344 1.552914 0.00015989949 -0.83034617 1.4173885 1.620925 0.00016350912 -0.81294478 1.4362241 1.555192 0.00016586488 -0.79514864 1.4511181 1.5562824 0.00016590867 -0.77696661 1.4677019 1.5249573 0.00016638401 -0.75843335 1.4815932 1.5680324 0.00016662967 -0.73948936 1.5027783 1.401677 0.00017062759 -0.72001765 1.5155065 1.4202473 0.00017421195 -0.7000796 1.5308187 1.3741596 0.00017593826" s = s + " -0.67965634 1.5432939 1.4005642 0.0001794182 -0.65874834 1.5628266 1.2723723 0.00018215045 -0.63728288 1.5767863 1.2949684 0.00018458866 -0.6152009 1.5943728 1.2109636 0.00018836483 -0.59249185 1.6124836 1.1000707 0.00019057883 -0.56903837 1.6211059 1.1553 0.00019433841 -0.54503602 1.6400816 1.1198099 0.00019489117 -0.52014209 1.6479358 1.2904679 0.00019787026 -0.49448261 1.6663249 1.23304 0.00020123336 -0.46782751 1.6726921 1.4163846 0.00020200546" s = s + " -0.44022652 1.6912618 1.2920816 0.00020616949 -0.41150914 1.7030679 1.3956081 0.00020791146 -0.3816396 1.7218323 1.2886276 0.00021143949 -0.35052328 1.7449335 1.2444879 0.00021608096 -0.31784963 1.7583988 1.3309349 0.00022137616 -0.28358669 1.7819671 1.1939449 0.0002282461 -0.24752936 1.8043432 1.1492554 0.0002334428 -0.20921434 1.8128613 1.4681372 0.00023912457 -0.16854843 1.8291524 1.6811082 0.00024296201 -0.1250677 1.8484229 1.8448104 0.00024956133" s = s + " -0.078452312 1.8935948 1.648603 0.00025539564 -0.02779895 1.9267952 1.7041255 0.00025829514 0.027706392 1.9619369 1.8949738 0.00026897317 0.089286662 1.9957209 2.2039115 0.00028266973 0.15856116 2.0391692 2.5297069 0.00030069415 0.23824409 2.0909362 2.9792799 0.00032333274 0.33270796 2.1330179 4.311973 0.00035243873 0.44958524 2.2800521 4.2114918 0.0003918326 0.60714202 2.457835 5.4652585 0.00044240428 0.64662561 2.5244955 5.5903544 0.00046150771" s = s + " 0.69039936 2.5833949 5.9261945 0.00048237641 0.73897923 2.6541195 6.3608976 0.00050414902 0.79395621 2.7338103 7.1352379 0.00054670677 0.85790385 2.8512874 8.0839936 0.00059870411 0.93404481 3.0496355 8.2089804 0.00066194575 1.0303788 3.2107748 10.398269 0.00075922721 1.1623571 3.4537569 14.242427 0.00090555619 1.3751885 4.0430952 19.400154 0.0011609585 1.5748041 4.717652 24.795177 0.001549166 1.8295984 5.9233078 29.935654 0.0023156824" s = s + " 2.0018937 6.5929082 40.856258 0.0030589721" v = strtoreal(tokens(s)) return(rowshape(v, 221)) } real matrix bu_urc_ct() { string scalar s real rowvector v s = "" s = s + " -5.1292249 -26.719405 -71.505058 -1813.2204 0.0033885902 -4.9772843 -23.023496 -106.2591 -824.40218 0.0026636323 -4.7675056 -18.937358 -106.17991 -284.37677 0.0017715683 -4.5953135 -16.764209 -68.301228 -462.16138 0.0012267776 -4.4175149 -14.124948 -60.038464 -323.3844 0.00093777691 -4.3074125 -12.945902 -44.683263 -336.47743 0.00079468331 -4.2275728 -12.030693 -36.566098 -352.38619 0.0006911358 -4.1638843 -11.350053 -32.774943 -292.7785 0.00063012787 -4.1111804 -10.738118 -32.793566 -222.42918 0.00058445069 -4.0655516 -10.355485 -24.941082 -269.16719 0.00055462829" s = s + " -4.025716 -9.9051615 -24.342994 -238.23165 0.00053346611 -3.9903382 -9.4980973 -24.469611 -202.11084 0.00051061829 -3.9579703 -9.1969574 -22.702701 -188.39198 0.00048941112 -3.8304297 -7.9968631 -15.371282 -186.22806 0.00041578111 -3.7359746 -7.1683122 -11.342941 -164.67362 0.00037278025 -3.6604348 -6.4961092 -9.8084411 -138.65435 0.00033483802 -3.5971452 -5.958473 -7.6698007 -139.80717 0.00031748485 -3.5420886 -5.5368953 -4.9209822 -151.27142 0.00029998883 -3.4933873 -5.161807 -3.8734324 -141.61612 0.00028230115 -3.4498312 -4.7945374 -4.5103958 -118.19984 0.00026593309" s = s + " -3.4098271 -4.5119686 -3.6630374 -109.96942 0.00026006414 -3.3731039 -4.253054 -2.6810539 -111.47716 0.00025015961 -3.3389269 -4.0161454 -2.2793812 -105.99232 0.00024376895 -3.3071958 -3.7837582 -2.6840288 -88.568817 0.00023631403 -3.2772698 -3.5989629 -1.4096929 -97.886557 0.00023105111 -3.2490638 -3.4004849 -1.8647137 -82.041248 0.00022374466 -3.2222496 -3.2298164 -1.6497471 -76.770146 0.00021599403 -3.1966701 -3.0914518 -0.46648607 -84.83216 0.00021510158 -3.1723541 -2.9442737 -0.25122362 -79.841118 0.00021196397 -3.1490499 -2.7984615 -0.60841107 -67.248271 0.00020392752" s = s + " -3.1266008 -2.6801993 0.18064468 -70.5668 0.00019867722 -3.1050614 -2.5556146 0.3048303 -66.331356 0.00019657521 -3.0842631 -2.4325536 0.055237751 -57.441312 0.00019363302 -3.064124 -2.3516058 1.9369931 -78.41441 0.00018938648 -3.0447284 -2.2387749 1.9473749 -76.524563 0.00018911144 -3.0258531 -2.138032 2.0895666 -75.536284 0.00018775206 -3.0075743 -2.0346668 1.9087771 -69.841588 0.00018362816 -2.9897721 -1.9406836 1.6726132 -61.191433 0.00018188085 -2.9724186 -1.8630055 2.1982439 -62.991422 0.00017953184 -2.9556152 -1.7764191 2.3943049 -63.323257 0.00017739996" s = s + " -2.939127 -1.7019814 2.9011393 -67.401446 0.00017631621 -2.9230915 -1.6194464 2.9100092 -65.734967 0.00017383308 -2.9073542 -1.5544466 3.3714972 -68.447922 0.00017132397 -2.8920421 -1.4746117 3.287292 -65.347757 0.00017065386 -2.8770072 -1.3999656 3.1762313 -62.357562 0.00016854506 -2.8622767 -1.3248669 2.9248375 -57.117555 0.00016775877 -2.8478144 -1.2576021 3.0275682 -57.490293 0.000168107 -2.8336368 -1.1832725 2.5774957 -50.417268 0.00016717041 -2.8196693 -1.1351672 3.2842425 -56.866936 0.00016507566 -2.8059229 -1.0897831 4.0987862 -65.845043 0.00016409938" s = s + " -2.7924244 -1.0357843 4.3139195 -66.451953 0.00016122474 -2.7791594 -0.98116276 4.4046139 -65.621184 0.00016021626 -2.7660918 -0.9356533 4.9513456 -71.803048 0.00016028087 -2.753295 -0.87603023 4.7595728 -68.431777 0.00016026079 -2.740577 -0.83901838 5.4578623 -74.961638 0.00015959372 -2.7281403 -0.7784162 4.9174947 -65.332097 0.00015821915 -2.7158486 -0.72739169 4.9094846 -64.184597 0.00015675119 -2.7037032 -0.67268668 4.4944629 -56.574353 0.00015596401 -2.691673 -0.62973528 4.5915487 -55.634219 0.00015596307 -2.6798653 -0.5769469 4.2675244 -50.038814 0.00015574682" s = s + " -2.6681478 -0.53657585 4.4389738 -50.652347 0.00015428007 -2.6566314 -0.48471837 4.0575016 -44.261354 0.00015527884 -2.6452172 -0.44006839 4.0839919 -44.101683 0.00015518621 -2.6339288 -0.3908084 3.592295 -35.487137 0.0001530364 -2.622692 -0.356808 3.9609952 -39.628037 0.00015324745 -2.6116268 -0.31585896 3.8918079 -37.337359 0.00015411009 -2.6007008 -0.27508044 3.8825003 -37.238795 0.00015260244 -2.5898604 -0.23595815 3.8374761 -35.703014 0.00015206038 -2.5791024 -0.19275104 3.6180677 -32.854471 0.00015233346 -2.568478 -0.14422428 2.9497493 -23.265165 0.00015187371" s = s + " -2.5578907 -0.11708115 3.4865837 -30.227654 0.00015146681 -2.5474479 -0.078928335 3.3787342 -28.511927 0.00015084094 -2.537051 -0.046002003 3.4911293 -29.808341 0.00015135005 -2.5267346 -0.0088433368 3.2582071 -25.573831 0.00015146718 -2.5165101 0.030239666 2.940012 -21.369329 0.00015071137 -2.5063486 0.057189111 3.1815506 -23.456797 0.00015069188 -2.4963173 0.096684467 2.8368008 -18.8104 0.00014896826 -2.4863552 0.13442678 2.5591894 -15.152288 0.00014803579 -2.4764134 0.16707005 2.3911364 -12.213498 0.00014642068 -2.4665145 0.19830071 2.2971835 -9.9591831 0.00014631912" s = s + " -2.4567076 0.22601828 2.4265859 -11.559708 0.00014551435 -2.4469298 0.25058604 2.6670623 -14.42383 0.00014500881 -2.4372232 0.27801061 2.6743891 -13.439708 0.0001446792 -2.4276224 0.31395462 2.2621507 -7.5122092 0.00014434048 -2.418023 0.3384695 2.3723621 -8.151206 0.00014468125 -2.408501 0.37196042 2.1008909 -5.0165956 0.00014364886 -2.3990018 0.40127084 1.850321 0.094670959 0.00014471542 -2.3895519 0.42357106 2.0564043 -2.0926551 0.00014384185 -2.3801178 0.44483141 2.3172042 -5.7502842 0.00014349328 -2.3707386 0.46799263 2.4510974 -7.8690664 0.00014289131" s = s + " -2.3613945 0.48453194 2.8496225 -12.430669 0.00014246123 -2.3521315 0.51224819 2.7172324 -10.805897 0.00014311951 -2.3428615 0.53499021 2.7393295 -10.47717 0.00014373853 -2.3336594 0.56658689 2.3561888 -5.6697405 0.0001437255 -2.324455 0.58537293 2.573145 -8.3767262 0.00014275271 -2.3152939 0.61050006 2.3544778 -4.5072384 0.00014221504 -2.3061545 0.63321402 2.2490796 -2.3168224 0.00014150262 -2.2970767 0.6571268 2.1628346 -1.1015245 0.00014192615 -2.2879866 0.67872982 2.0570538 1.1234809 0.00014112228 -2.2789051 0.68916911 2.6927433 -7.4046893 0.00014024701" s = s + " -2.2698488 0.70828441 2.7444158 -7.4737131 0.00013996383 -2.2608368 0.73335832 2.4663581 -3.2373265 0.00014066022 -2.2518533 0.75221949 2.5753942 -4.6643976 0.00014096756 -2.2428768 0.77080224 2.7153251 -6.6507102 0.00013965135 -2.233945 0.79382399 2.6031705 -5.2426895 0.00013841971 -2.2250064 0.81629951 2.4828676 -3.8413947 0.00013786648 -2.2160811 0.84061753 2.2139888 -0.016662817 0.00013701353 -2.2071578 0.85684905 2.4074056 -2.462436 0.00013804943 -2.19826 0.87649237 2.4263314 -3.0569545 0.00013797874 -2.1893693 0.89717147 2.3049354 -0.91761864 0.00013818453" s = s + " -2.1804723 0.9140739 2.3926728 -2.1744109 0.00013840677 -2.1715633 0.93349657 2.3061015 -0.90584152 0.00013761597 -2.1626276 0.93971963 2.8345674 -7.1659328 0.00013832364 -2.1537543 0.95788791 2.9207276 -8.9025672 0.00013807326 -2.1448725 0.97685127 2.8226055 -7.3310995 0.00013747915 -2.1360234 1.0018044 2.5081261 -3.7834633 0.00013608001 -2.127141 1.0221381 2.3591607 -1.755275 0.00013647218 -2.1182572 1.0376907 2.4550352 -2.8045144 0.00013606637 -2.1093539 1.0545503 2.4832597 -3.6485319 0.00013663996 -2.1004392 1.0707521 2.516962 -4.3821238 0.00013792304" s = s + " -2.0915177 1.076621 3.1197524 -12.823056 0.00013844332 -2.0825577 1.0893451 3.1880812 -13.158546 0.00013863815 -2.0736342 1.1094116 2.96821 -10.46919 0.00013787602 -2.0646655 1.1188882 3.388107 -16.443775 0.00013929772 -2.0556978 1.1354394 3.3064686 -15.043946 0.00013896461 -2.0467186 1.154045 3.1419242 -12.67585 0.00014019431 -2.0377119 1.1701541 3.1345858 -12.88273 0.0001403696 -2.0286865 1.1840489 3.2405907 -14.277327 0.00014100853 -2.0196066 1.1908254 3.6997872 -20.796341 0.0001399728 -2.0105035 1.1986852 4.0948873 -26.41497 0.0001400136" s = s + " -2.0014101 1.2174824 3.838716 -22.856574 0.000139146 -1.9922608 1.225849 4.1415715 -26.562092 0.00013826548 -1.9831055 1.242126 4.1281558 -26.661964 0.00013867907 -1.9739381 1.2638568 3.7981407 -22.613934 0.00013827443 -1.9646817 1.2662124 4.4468385 -31.496535 0.00013817797 -1.9554013 1.2818893 4.2857828 -28.94869 0.00013691489 -1.9460905 1.2949651 4.3324966 -29.752916 0.00013695602 -1.9367325 1.30594 4.5148127 -32.156138 0.0001365365 -1.9273449 1.3196473 4.4970527 -32.103802 0.00013685557 -1.917942 1.3390194 4.2389864 -28.474253 0.0001368387" s = s + " -1.9084974 1.3616852 3.8076995 -23.024484 0.00013672194 -1.8989633 1.3792449 3.5856095 -19.966922 0.00013753257 -1.8893909 1.4027174 3.0218952 -12.727231 0.00013828068 -1.8797154 1.4114457 3.2282461 -14.933104 0.00013922397 -1.8700179 1.4268665 3.1472097 -14.232556 0.00014006429 -1.8602368 1.4356002 3.4593436 -18.790799 0.00014016004 -1.8504258 1.4500848 3.3863764 -17.434797 0.00014059273 -1.8405255 1.4662367 3.2540039 -16.131579 0.0001412837 -1.8305438 1.4807586 3.1937663 -15.461995 0.00014214132 -1.8204991 1.4948178 3.2497226 -16.927232 0.00014075805" s = s + " -1.810351 1.5097792 3.1988594 -16.425651 0.00014139416 -1.8001244 1.5250088 3.065401 -14.728335 0.00014027316 -1.789788 1.5308619 3.408534 -18.667621 0.00014037122 -1.7793754 1.5397714 3.6001242 -20.58658 0.00014081775 -1.7688669 1.5515074 3.6250524 -20.123236 0.00014178541 -1.7582836 1.5657976 3.710857 -21.530377 0.00014184305 -1.7476215 1.5888913 3.3709576 -18.128033 0.00014275169 -1.7368102 1.6075322 3.0628266 -13.698175 0.00014441364 -1.7258802 1.6217158 2.9465403 -10.831516 0.00014414746 -1.7147789 1.6263007 3.3347902 -14.567624 0.00014487648" s = s + " -1.7036055 1.6385091 3.5964522 -19.478659 0.0001458091 -1.6923035 1.6575011 3.387545 -16.471074 0.00014675254 -1.6808178 1.6735087 3.3470112 -15.979085 0.00014779663 -1.6692208 1.6976318 2.9191207 -11.105685 0.00014845919 -1.6574499 1.712672 3.0190391 -12.544351 0.00014911957 -1.6455203 1.7321196 2.7843557 -8.8951804 0.00015133137 -1.6334202 1.7505771 2.7140815 -8.4134683 0.0001511422 -1.6211175 1.7647263 2.830155 -9.3961723 0.00015302209 -1.6086311 1.7808217 2.8613978 -9.3538199 0.00015394567 -1.5959047 1.7966812 2.9838201 -11.103841 0.00015518715" s = s + " -1.5830072 1.8192877 2.6503963 -5.56307 0.00015500475 -1.5698372 1.8341156 2.739595 -5.1385562 0.00015687623 -1.5564282 1.8562127 2.501214 -0.78716052 0.00015882836 -1.5427607 1.8773714 2.3435865 2.8472436 0.00016039372 -1.5288002 1.8956652 2.2980392 5.6595326 0.00016232841 -1.5145361 1.9093733 2.6414426 2.3517865 0.00016365345 -1.4999657 1.9260694 2.8047512 1.689443 0.00016566139 -1.4850839 1.9527455 2.5169865 7.6460826 0.00016641569 -1.4698449 1.9750883 2.5025928 9.8837033 0.00017034911 -1.4541516 1.9847568 3.1504806 4.2753172 0.00017141782" s = s + " -1.438104 2.0099117 3.2157148 4.8328179 0.00017448019 -1.4215628 2.0239976 3.8039014 -0.28153304 0.00017656234 -1.4046128 2.0543491 3.5563272 5.7751371 0.0001769942 -1.38714 2.0916555 3.0732965 14.767552 0.00017892036 -1.3690755 2.1113911 3.4866191 12.462613 0.00018268588 -1.3503417 2.1236448 4.3946426 2.4200196 0.00018601054 -1.3310719 2.1629672 4.0026267 9.9316285 0.000187214 -1.3111432 2.1994885 4.0790645 10.240756 0.00019061032 -1.2904021 2.2420436 3.7791777 15.485402 0.00019388647 -1.2689099 2.3008585 3.000574 25.4151 0.00019911109" s = s + " -1.2464024 2.32862 3.7709529 17.058953 0.00020395617 -1.2228862 2.3634295 4.2775097 11.250807 0.0002094463 -1.1982425 2.4090033 4.3932276 8.9032503 0.00021405604 -1.1723118 2.4458887 4.756895 6.8177565 0.00022029962 -1.1450172 2.4912967 5.0152698 2.9440321 0.00022563475 -1.1161953 2.5419524 5.0310148 2.6089919 0.00023322639 -1.0856992 2.6004777 4.6951057 7.8180145 0.00024328614 -1.0532375 2.684496 3.174349 24.984198 0.00025189864 -1.0185068 2.760494 2.2321465 34.995579 0.00026108219 -0.98105983 2.8205207 1.8424027 42.493438 0.00027192017" s = s + " -0.94030112 2.8636585 2.365293 34.26805 0.00028174168 -0.89611941 2.9410406 1.5984543 44.446099 0.00029244718 -0.84734516 3.0284974 0.36655117 60.86628 0.00030364579 -0.79301687 3.1032474 -0.070383002 67.919785 0.0003225959 -0.73134643 3.1699769 0.57238705 56.455304 0.00034359911 -0.65988222 3.2397533 1.0038603 55.590683 0.00036281585 -0.57498908 3.3704549 -0.86853477 90.162869 0.00040371043 -0.46882616 3.4565692 1.0856959 73.944809 0.00045839375 -0.32542503 3.6467386 0.97278076 80.096055 0.00053559421 -0.28927049 3.6667954 2.6798444 58.498617 0.00054654176" s = s + " -0.24936161 3.7425553 1.9014721 69.517502 0.00056565084 -0.20449854 3.8183521 1.2114951 81.493244 0.00059004685 -0.15389825 3.8686285 3.7704385 48.132614 0.00063012289 -0.094616438 3.9007553 8.2410287 -17.45936 0.00066143145 -0.023862549 4.0207808 7.1436152 20.57007 0.00075075038 0.065520718 4.1594945 6.7479375 46.321645 0.00082719858 0.1883404 4.3459788 9.7447054 8.0854981 0.00096268852 0.38863151 4.6531364 14.968841 -8.3553453 0.0013263463 0.57587025 5.070594 21.732595 -63.042215 0.0017499956 0.81904547 5.2559469 56.007846 -411.78944 0.0026299291" s = s + " 0.97720956 6.6104917 32.425153 -157.35229 0.0034296861" v = strtoreal(tokens(s)) return(rowshape(v, 221)) } real matrix bu_urc_probs() { string scalar s real rowvector v s = "" s = s + " 1e-04 -3.71901649 2e-04 -3.5400838 5e-04 -3.29052673 0.001 -3.09023231 0.002 -2.87816174 0.003 -2.74778139 0.004 -2.65206981 0.005 -2.5758293 0.006 -2.51214433 0.007 -2.45726339" s = s + " 0.008 -2.40891555 0.009 -2.36561813 0.01 -2.32634787 0.015 -2.17009038 0.02 -2.05374891 0.025 -1.95996398 0.03 -1.88079361 0.035 -1.81191067 0.04 -1.75068607 0.045 -1.69539771" s = s + " 0.05 -1.64485363 0.055 -1.59819314 0.06 -1.55477359 0.065 -1.51410189 0.07 -1.47579103 0.075 -1.43953147 0.08 -1.40507156 0.085 -1.37220381 0.09 -1.34075503 0.095 -1.31057911" s = s + " 0.1 -1.28155157 0.105 -1.25356544 0.11 -1.22652812 0.115 -1.20035886 0.12 -1.17498679 0.125 -1.15034938 0.13 -1.12639113 0.135 -1.10306256 0.14 -1.08031934 0.145 -1.05812162" s = s + " 0.15 -1.03643339 0.155 -1.01522203 0.16 -0.99445788 0.165 -0.97411388 0.17 -0.95416525 0.175 -0.93458929 0.18 -0.91536509 0.185 -0.89647336 0.19 -0.8778963 0.195 -0.85961736" s = s + " 0.2 -0.84162123 0.205 -0.82389363 0.21 -0.80642125 0.215 -0.78919165 0.22 -0.77219321 0.225 -0.75541503 0.23 -0.73884685 0.235 -0.72247905 0.24 -0.70630256 0.245 -0.69030882" s = s + " 0.25 -0.67448975 0.255 -0.65883769 0.26 -0.64334541 0.265 -0.62800601 0.27 -0.61281299 0.275 -0.59776013 0.28 -0.58284151 0.285 -0.5680515 0.29 -0.55338472 0.295 -0.53883603" s = s + " 0.3 -0.52440051 0.305 -0.51007346 0.31 -0.49585035 0.315 -0.48172685 0.32 -0.4676988 0.325 -0.45376219 0.33 -0.43991317 0.335 -0.42614801 0.34 -0.41246313 0.345 -0.39885507" s = s + " 0.35 -0.38532047 0.355 -0.37185609 0.36 -0.35845879 0.365 -0.34512553 0.37 -0.33185335 0.375 -0.31863936 0.38 -0.30548079 0.385 -0.2923749 0.39 -0.27931903 0.395 -0.26631061" s = s + " 0.4 -0.2533471 0.405 -0.24042603 0.41 -0.22754498 0.415 -0.21470157 0.42 -0.20189348 0.425 -0.18911843 0.43 -0.17637416 0.435 -0.16365849 0.44 -0.15096922 0.445 -0.13830421" s = s + " 0.45 -0.12566135 0.455 -0.11303854 0.46 -0.10043372 0.465 -0.08784484 0.47 -0.07526986 0.475 -0.06270678 0.48 -0.05015358 0.485 -0.03760829 0.49 -0.02506891 0.495 -0.01253347" s = s + " 0.5 0 0.505 0.01253347 0.51 0.02506891 0.515 0.03760829 0.52 0.05015358 0.525 0.06270678 0.53 0.07526986 0.535 0.08784484 0.54 0.10043372 0.545 0.11303854" s = s + " 0.55 0.12566135 0.555 0.13830421 0.56 0.15096922 0.565 0.16365849 0.57 0.17637416 0.575 0.18911843 0.58 0.20189348 0.585 0.21470157 0.59 0.22754498 0.595 0.24042603" s = s + " 0.6 0.2533471 0.605 0.26631061 0.61 0.27931903 0.615 0.2923749 0.62 0.30548079 0.625 0.31863936 0.63 0.33185335 0.635 0.34512553 0.64 0.35845879 0.645 0.37185609" s = s + " 0.65 0.38532047 0.655 0.39885507 0.66 0.41246313 0.665 0.42614801 0.67 0.43991317 0.675 0.45376219 0.68 0.4676988 0.685 0.48172685 0.69 0.49585035 0.695 0.51007346" s = s + " 0.7 0.52440051 0.705 0.53883603 0.71 0.55338472 0.715 0.5680515 0.72 0.58284151 0.725 0.59776013 0.73 0.61281299 0.735 0.62800601 0.74 0.64334541 0.745 0.65883769" s = s + " 0.75 0.67448975 0.755 0.69030882 0.76 0.70630256 0.765 0.72247905 0.77 0.73884685 0.775 0.75541503 0.78 0.77219321 0.785 0.78919165 0.79 0.80642125 0.795 0.82389363" s = s + " 0.8 0.84162123 0.805 0.85961736 0.81 0.8778963 0.815 0.89647336 0.82 0.91536509 0.825 0.93458929 0.83 0.95416525 0.835 0.97411388 0.84 0.99445788 0.845 1.01522203" s = s + " 0.85 1.03643339 0.855 1.05812162 0.86 1.08031934 0.865 1.10306256 0.87 1.12639113 0.875 1.15034938 0.88 1.17498679 0.885 1.20035886 0.89 1.22652812 0.895 1.25356544" s = s + " 0.9 1.28155157 0.905 1.31057911 0.91 1.34075503 0.915 1.37220381 0.92 1.40507156 0.925 1.43953147 0.93 1.47579103 0.935 1.51410189 0.94 1.55477359 0.945 1.59819314" s = s + " 0.95 1.64485363 0.955 1.69539771 0.96 1.75068607 0.965 1.81191067 0.97 1.88079361 0.975 1.95996398 0.98 2.05374891 0.985 2.17009038 0.99 2.32634787 0.991 2.36561813" s = s + " 0.992 2.40891555 0.993 2.45726339 0.994 2.51214433 0.995 2.5758293 0.996 2.65206981 0.997 2.74778139 0.998 2.87816174 0.999 3.09023231 0.9995 3.29052673 0.9998 3.5400838" s = s + " 0.9999 3.71901649" v = strtoreal(tokens(s)) return(rowshape(v, 221)) } // MacKinnon (1996) p-value for a unit-root t-statistic. // cs: "nc","c","ct" ; nobs: sample size for the response surface (0 = asymptotic) real scalar bu_fpval(real scalar stat, string scalar cs, real scalar nobs) { real matrix beta, pc, X, W real colvector crits, probs, cnorm, wght, yv, g real scalar i, np, nph, ic, np1, np2, nvar, eta, z if (cs=="nc") { beta = bu_urc_nc() nvar = 3 } else if (cs=="c") { beta = bu_urc_c() nvar = 3 } else { beta = bu_urc_ct() nvar = 4 } pc = bu_urc_probs() probs = pc[.,1] cnorm = pc[.,2] wght = beta[., cols(beta)] crits = beta[.,1] if (nobs > 0) { eta = 1/nobs crits = crits + beta[.,2]*eta + beta[.,3]*(eta*eta) if (nvar==4) crits = crits + beta[.,4]*(eta*eta*eta) } np = 9 nph = 4 ic = 222 for (i=1; i<=221; i++) { if (crits[i] > stat) { ic = i i = 222 } } np1 = ic - nph if (np1 < 1) np1 = 1 np2 = np1 + np - 1 if (np2 > 221) { np2 = 221 np1 = np2 - np + 1 } X = J(np,4,1) yv = J(np,1,0) W = J(np,1,0) for (i=1; i<=np; i++) { X[i,2] = crits[np1+i-1] X[i,3] = crits[np1+i-1]*crits[np1+i-1] X[i,4] = crits[np1+i-1]*crits[np1+i-1]*crits[np1+i-1] yv[i] = cnorm[np1+i-1] W[i] = 1/(wght[np1+i-1]*wght[np1+i-1]) } g = invsym(cross(X,W,X)) * cross(X,W,yv) z = g[1] + g[2]*stat + g[3]*stat*stat + g[4]*stat*stat*stat return(normal(z)) } // <<>> end