*! lsmlm 1.1.0 20jul2026 Ibrahim Ongoren, Pamukkale University *! Lee-Strazicich minimum LM unit root test with zero, one, or two breaks. *! Models: crash (level shift) and break (level + trend shift). breaks(0) is *! the Schmidt-Phillips no-break LM test. *! *! General-purpose command: it operates on the current tsset time series *! (or a series given via timevar()) and can be restricted to any sub-sample *! with if / in. The test statistic is computed in native Stata (adapted from *! the author's own ls1break_lm / ls2break_lm routines). *! *! Critical values: the interpolated CV tables (interpolation over the sample *! size T and, for the trend-break model, over the break fraction lambda) are *! ported from the -leestra- package by H. Ozan Eruygur (SSC s459688), which *! is distributed under the MIT License. The tables ultimately come from *! Schmidt & Phillips (1992) and Lee & Strazicich (2003, 2013). * * Method (Lee & Strazicich 2003, 2013; Schmidt & Phillips 1992): * 1) For each candidate break date (pair), regress D.y on the differenced * deterministic terms DZ to recover the shift magnitudes, then build the * LM-detrended series S(t) (S(1)=0, D.S = residual of D.y on DZ). * 2) Regress D.y on DZ, S(t-1) and lagged D.S terms; general-to-specific * lag selection from maxlag() downward (keep the highest lag with * |t| >= 1.645). * 3) Pick the break date(s) that MINIMISE the t-statistic on S(t-1). * 4) Compare that minimum statistic to the interpolated LS/SP critical values. * * Deterministic terms by model (per break j at index b): * impulse B_j = 1{t == b+1}, step D_j = 1{t > b}, trend T_j = (t-b)*1{t>b} * crash model : DZ regressors = {const, B_j}; S subtracts level steps * break model : DZ regressors = {const, B_j, D_j}; S subtracts level + trend * no break : DZ regressors = {const} program define lsmlm, rclass version 14.0 syntax varname(numeric) [if] [in] , /// [ Breaks(integer 2) /// Model(string) /// MAXLag(integer -1) /// TRim(real 0.10) /// MINDist(integer -1) /// TIMEvar(varname numeric) /// noPRINT ] * ---------- validate options ---------- if "`model'" == "" local model "break" if !inlist("`model'", "break", "crash") { di as err "model() must be break or crash" exit 198 } local mflag = cond("`model'" == "crash", 1, 2) if !inlist(`breaks', 0, 1, 2) { di as err "breaks() must be 0, 1, or 2" exit 198 } if `trim' <= 0 | `trim' >= 0.5 { di as err "trim() must be in (0, 0.5)" exit 198 } * ---------- determine the time variable ---------- if "`timevar'" == "" { cap tsset if _rc { di as err "data must be tsset, or specify timevar(), before using lsmlm" exit 198 } local timevar `r(timevar)' } marksample touse qui replace `touse' = 0 if missing(`varlist') | missing(`timevar') preserve qui keep if `touse' sort `timevar' cap qui isid `timevar' if _rc { di as err "the time variable is not unique in the current sample." di as err "lsmlm operates on a single time series: restrict to one panel" di as err "unit with if/in (e.g. if id==1), or tsset a single series." exit 198 } qui count local T = r(N) if `T' < 20 { di as err "not enough observations after if/in and missing (need ~20+)" exit 198 } if `maxlag' < 0 { local maxlag = floor(4*(`T'/100)^0.25) if `maxlag' < 1 local maxlag = 1 } if `mindist' < 0 local mindist = 2 tempvar yy tt dy qui gen double `yy' = `varlist' qui gen long `tt' = _n qui gen double `dy' = `yy' - `yy'[_n-1] local lo = ceil(`trim' * `T') local hi = floor((1 - `trim') * `T') if `lo' < 3 local lo = 3 if `hi' > `T' - 3 local hi = `T' - 3 * Collinearity guard: keep the earliest candidate at least maxlag+2 so a * level dummy cannot be constant across the lagged-regression sample and * collide with the constant. if `lo' < `maxlag' + 2 local lo = `maxlag' + 2 if `breaks' > 0 & `lo' > `hi' { di as err "sample too short for the chosen trim() and maxlag(): empty break window" exit 198 } * ---------- best-configuration storage ---------- local best_tau = . local best_k = . local best_b1 = . local best_b2 = . local best_year1 = . local best_year2 = . local best_lam1 = . local best_lam2 = . local best_tB1 = . local best_tB2 = . local best_tD1 = . local best_tD2 = . local best_sig1 = . local best_sig2 = . local best_N = . if `breaks' == 0 { * ===================== NO BREAK (Schmidt-Phillips) ===================== capture drop _ls_S _ls_dS _ls_LS capture drop _ls_LdS* capture quietly regress `dy' if `tt' > 1 if !_rc { local btrend = _b[_cons] local y1 = `yy'[1] local psi = `y1' - (`btrend' * 1) qui gen double _ls_S = `yy' - `psi' - (`btrend' * `tt') qui gen double _ls_dS = _ls_S - _ls_S[_n-1] qui gen double _ls_LS = _ls_S[_n-1] local done = 0 forvalues kk = `maxlag'(-1)0 { capture drop _ls_LdS* local laglist "" if `kk' > 0 { forvalues L = 1/`kk' { qui gen double _ls_LdS`L' = _ls_dS[_n-`L'] local laglist "`laglist' _ls_LdS`L'" } } capture quietly regress `dy' _ls_LS `laglist' if !_rc & e(N) > 20 { local tau = _b[_ls_LS] / _se[_ls_LS] if `kk' == 0 { local best_k = 0 local best_tau = `tau' local best_N = e(N) local done = 1 continue, break } else { local lastvar "_ls_LdS`kk'" local tlast = abs(_b[`lastvar'] / _se[`lastvar']) if `tlast' >= 1.645 { local best_k = `kk' local best_tau = `tau' local best_N = e(N) local done = 1 continue, break } } } } } } else if `breaks' == 1 { * ===================== ONE BREAK ===================== forvalues b1 = `lo'/`hi' { capture drop _ls_D1 _ls_T1 _ls_B1 _ls_S _ls_dS _ls_LS capture drop _ls_LdS* qui gen double _ls_D1 = (`tt' > `b1') qui gen double _ls_T1 = (`tt' - `b1') * (`tt' > `b1') qui gen double _ls_B1 = (`tt' == `b1' + 1) if `mflag' == 2 local dvars "_ls_B1 _ls_D1" else local dvars "_ls_B1" capture quietly regress `dy' `dvars' if `tt' > 1 if !_rc { local btrend = _b[_cons] local blev1 = _b[_ls_B1] local y1 = `yy'[1] local psi = `y1' - (`btrend' * 1) if `mflag' == 2 { local btr1 = _b[_ls_D1] qui gen double _ls_S = `yy' - `psi' - (`btrend' * `tt') - (`blev1' * _ls_D1) - (`btr1' * _ls_T1) } else { qui gen double _ls_S = `yy' - `psi' - (`btrend' * `tt') - (`blev1' * _ls_D1) } qui gen double _ls_dS = _ls_S - _ls_S[_n-1] qui gen double _ls_LS = _ls_S[_n-1] local final_k = . local final_tau = . local final_N = . local final_tB1 = . local final_tD1 = . local done = 0 forvalues kk = `maxlag'(-1)0 { capture drop _ls_LdS* local laglist "" if `kk' > 0 { forvalues L = 1/`kk' { qui gen double _ls_LdS`L' = _ls_dS[_n-`L'] local laglist "`laglist' _ls_LdS`L'" } } capture quietly regress `dy' `dvars' _ls_LS `laglist' if !_rc & e(N) > 20 { local tau = _b[_ls_LS] / _se[_ls_LS] local ok = 0 if `kk' == 0 local ok = 1 else { local lastvar "_ls_LdS`kk'" if abs(_b[`lastvar'] / _se[`lastvar']) >= 1.645 local ok = 1 } if `ok' { local final_k = `kk' local final_tau = `tau' local final_N = e(N) local final_tB1 = _b[_ls_B1] / _se[_ls_B1] if `mflag' == 2 local final_tD1 = _b[_ls_D1] / _se[_ls_D1] local done = 1 continue, break } } } if `done' == 1 { if missing(`best_tau') | `final_tau' < `best_tau' { local best_tau = `final_tau' local best_k = `final_k' local best_b1 = `b1' local best_year1 = `timevar'[`b1'] local best_lam1 = `b1' / `T' local best_tB1 = `final_tB1' local best_tD1 = `final_tD1' if `mflag' == 2 local best_sig1 = (max(abs(`final_tB1'), abs(`final_tD1')) >= 1.645) else local best_sig1 = (abs(`final_tB1') >= 1.645) local best_N = `final_N' } } } capture drop _ls_D1 _ls_T1 _ls_B1 _ls_S _ls_dS _ls_LS capture drop _ls_LdS* } } else { * ===================== TWO BREAKS ===================== forvalues b1 = `lo'/`hi' { forvalues b2 = `lo'/`hi' { if (`b2' > `b1' + `mindist') { capture drop _ls_D1 _ls_D2 _ls_T1 _ls_T2 _ls_B1 _ls_B2 _ls_S _ls_dS _ls_LS capture drop _ls_LdS* qui gen double _ls_D1 = (`tt' > `b1') qui gen double _ls_D2 = (`tt' > `b2') qui gen double _ls_T1 = (`tt' - `b1') * (`tt' > `b1') qui gen double _ls_T2 = (`tt' - `b2') * (`tt' > `b2') qui gen double _ls_B1 = (`tt' == `b1' + 1) qui gen double _ls_B2 = (`tt' == `b2' + 1) if `mflag' == 2 local dvars "_ls_B1 _ls_B2 _ls_D1 _ls_D2" else local dvars "_ls_B1 _ls_B2" capture quietly regress `dy' `dvars' if `tt' > 1 if !_rc { local btrend = _b[_cons] local blev1 = _b[_ls_B1] local blev2 = _b[_ls_B2] local y1 = `yy'[1] local psi = `y1' - (`btrend' * 1) if `mflag' == 2 { local btr1 = _b[_ls_D1] local btr2 = _b[_ls_D2] qui gen double _ls_S = `yy' - `psi' - (`btrend' * `tt') - (`blev1' * _ls_D1) - (`blev2' * _ls_D2) - (`btr1' * _ls_T1) - (`btr2' * _ls_T2) } else { qui gen double _ls_S = `yy' - `psi' - (`btrend' * `tt') - (`blev1' * _ls_D1) - (`blev2' * _ls_D2) } qui gen double _ls_dS = _ls_S - _ls_S[_n-1] qui gen double _ls_LS = _ls_S[_n-1] local final_k = . local final_tau = . local final_N = . local final_tB1 = . local final_tB2 = . local final_tD1 = . local final_tD2 = . local done = 0 forvalues kk = `maxlag'(-1)0 { capture drop _ls_LdS* local laglist "" if `kk' > 0 { forvalues L = 1/`kk' { qui gen double _ls_LdS`L' = _ls_dS[_n-`L'] local laglist "`laglist' _ls_LdS`L'" } } capture quietly regress `dy' `dvars' _ls_LS `laglist' if !_rc & e(N) > 20 { local tau = _b[_ls_LS] / _se[_ls_LS] local ok = 0 if `kk' == 0 local ok = 1 else { local lastvar "_ls_LdS`kk'" if abs(_b[`lastvar'] / _se[`lastvar']) >= 1.645 local ok = 1 } if `ok' { local final_k = `kk' local final_tau = `tau' local final_N = e(N) local final_tB1 = _b[_ls_B1] / _se[_ls_B1] local final_tB2 = _b[_ls_B2] / _se[_ls_B2] if `mflag' == 2 { local final_tD1 = _b[_ls_D1] / _se[_ls_D1] local final_tD2 = _b[_ls_D2] / _se[_ls_D2] } local done = 1 continue, break } } } if `done' == 1 { if missing(`best_tau') | `final_tau' < `best_tau' { local best_tau = `final_tau' local best_k = `final_k' local best_b1 = `b1' local best_b2 = `b2' local best_year1 = `timevar'[`b1'] local best_year2 = `timevar'[`b2'] local best_lam1 = `b1' / `T' local best_lam2 = `b2' / `T' local best_tB1 = `final_tB1' local best_tB2 = `final_tB2' local best_tD1 = `final_tD1' local best_tD2 = `final_tD2' if `mflag' == 2 { local best_sig1 = (max(abs(`final_tB1'), abs(`final_tD1')) >= 1.645) local best_sig2 = (max(abs(`final_tB2'), abs(`final_tD2')) >= 1.645) } else { local best_sig1 = (abs(`final_tB1') >= 1.645) local best_sig2 = (abs(`final_tB2') >= 1.645) } local best_N = `final_N' } } } capture drop _ls_D1 _ls_D2 _ls_T1 _ls_T2 _ls_B1 _ls_B2 _ls_S _ls_dS _ls_LS capture drop _ls_LdS* } } } } * ---------- did we find an estimable configuration? ---------- if missing(`best_tau') { di as err "no estimable configuration was found; try a larger sample," di as err "a smaller maxlag(), or a smaller trim()." exit 498 } * ---------- break fraction(s): leestra's lag-adjusted convention ---------- * lambda = (break_index - (k+2)) / (T - (k+2) + 1), matching -leestra- so the * CV lookup lands on the same table cell as the reference implementation. if `breaks' >= 1 { local lower_eff = `best_k' + 2 local best_lam1 = (`best_b1' - `lower_eff') / (`T' - `lower_eff' + 1) if `breaks' == 2 { local best_lam2 = (`best_b2' - `lower_eff') / (`T' - `lower_eff' + 1) } } * ---------- interpolated critical values (Mata) ---------- local lam1arg = cond(missing(`best_lam1'), 0, `best_lam1') local lam2arg = cond(missing(`best_lam2'), 0, `best_lam2') mata: lsmlm_cvlookup(`breaks', `mflag', `lam1arg', `lam2arg', `best_N') local cv1 = r(cv1) local cv5 = r(cv5) local cv10 = r(cv10) local reject10 = (`best_tau' <= `cv10') local reject5 = (`best_tau' <= `cv5') local reject1 = (`best_tau' <= `cv1') local siglevel "Not significant" if `reject10' local siglevel "10%" if `reject5' local siglevel "5%" if `reject1' local siglevel "1%" * sample range in time units (before dropping the working copy) local tmin = `timevar'[1] local tmax = `timevar'[`T'] local fmt : format `timevar' if "`fmt'" == "" local fmt "%9.0g" local tmintxt = string(`tmin', "`fmt'") local tmaxtxt = string(`tmax', "`fmt'") if `breaks' >= 1 local b1txt = string(`best_year1', "`fmt'") if `breaks' == 2 local b2txt = string(`best_year2', "`fmt'") restore * ---------- returned results ---------- return scalar tau = `best_tau' return scalar k = `best_k' return scalar breaks = `breaks' return scalar N = `best_N' return scalar maxlag = `maxlag' return scalar trim = `trim' return scalar cv1 = `cv1' return scalar cv5 = `cv5' return scalar cv10 = `cv10' return scalar reject1 = `reject1' return scalar reject5 = `reject5' return scalar reject10 = `reject10' if `breaks' >= 1 { return scalar tb1 = `best_year1' return scalar tb1_index = `best_b1' return scalar lambda1 = `best_lam1' return scalar break1_sig10 = `best_sig1' return scalar tB1 = `best_tB1' if `mflag' == 2 return scalar tD1 = `best_tD1' } if `breaks' == 2 { return scalar mindist = `mindist' return scalar tb2 = `best_year2' return scalar tb2_index = `best_b2' return scalar lambda2 = `best_lam2' return scalar break2_sig10 = `best_sig2' return scalar tB2 = `best_tB2' if `mflag' == 2 return scalar tD2 = `best_tD2' } return local siglevel "`siglevel'" return local model "`model'" return local timevar "`timevar'" return local varname "`varlist'" * ---------- display ---------- if "`print'" == "" { if `reject10' local decision "Reject unit root at `siglevel'" else local decision "Do not reject unit root" local mdesc = cond(`mflag'==2, "trend break (level + trend shift)", "crash (level shift)") if `breaks' == 0 local mdesc "no break (Schmidt-Phillips)" di di as txt "{hline 66}" di as res " Lee-Strazicich minimum LM unit root test" di as txt "{hline 66}" di as txt " Series: " as res "`varlist'" di as txt " Model: " as res "`mdesc'" di as txt " Breaks: " as res "`breaks'" as txt " (endogenously selected)" di as txt " Max lag: " as res "`maxlag'" as txt " Lags chosen (GTOS): " as res "`best_k'" di as txt " Obs (regr.): " as res "`best_N'" di as txt " Sample: " as res "`tmintxt'" as txt " to " as res "`tmaxtxt'" if `breaks' >= 1 { di as txt " Break 1: " as res "`b1txt'" as txt " (lambda1 = " as res %5.3f `best_lam1' as txt ")" } if `breaks' == 2 { di as txt " Break 2: " as res "`b2txt'" as txt " (lambda2 = " as res %5.3f `best_lam2' as txt ")" } di as txt "{hline 66}" di as txt " Test statistic (t on S(t-1)): " as res %9.4f `best_tau' di as txt " Critical values: 1% " as res %7.3f `cv1' as txt " 5% " as res %7.3f `cv5' as txt " 10% " as res %7.3f `cv10' di as txt " Decision: " as res "`decision'" di as txt "{hline 66}" if `breaks' >= 1 { di as txt " Break significance (|t| on dummies >= 1.645):" local s1txt = cond(`best_sig1', "significant", "not significant") di as txt " Break 1: " as res "`s1txt'" if `breaks' == 2 { local s2txt = cond(`best_sig2', "significant", "not significant") di as txt " Break 2: " as res "`s2txt'" } di as txt "{hline 66}" } } end *=========================================================================== * Mata: interpolated critical-value engine. * * The CV tables and interpolation routines below are ported from the -leestra- * package by H. Ozan Eruygur (SSC s459688), distributed under the MIT License. * Underlying sources: Schmidt & Phillips (1992); Lee & Strazicich (2003, 2013). *=========================================================================== version 14.0 mata: mata set matastrict off real rowvector lsmlm_nlookup() { return((100, 250, 500, 1000)) } // No break (Schmidt-Phillips) real matrix lsmlm_cv0() { real matrix M M = J(4, 3, .) M[1,1] = -3.597; M[1,2] = -3.031; M[1,3] = -2.745 M[2,1] = -3.572; M[2,2] = -3.023; M[2,3] = -2.747 M[3,1] = -3.570; M[3,2] = -3.021; M[3,3] = -2.748 M[4,1] = -3.566; M[4,2] = -3.023; M[4,3] = -2.748 return(M) } // One break, crash model real matrix lsmlm_cvc1() { real matrix M M = J(4, 3, .) M[1,1] = -4.084; M[1,2] = -3.487; M[1,3] = -3.185 M[2,1] = -3.987; M[2,2] = -3.387; M[2,3] = -3.076 M[3,1] = -3.840; M[3,2] = -3.277; M[3,3] = -2.985 M[4,1] = -3.798; M[4,2] = -3.230; M[4,3] = -2.925 return(M) } // Two breaks, crash model real matrix lsmlm_cvc2() { real matrix M M = J(4, 3, .) M[1,1] = -4.073; M[1,2] = -3.563; M[1,3] = -3.296 M[2,1] = -4.101; M[2,2] = -3.594; M[2,3] = -3.345 M[3,1] = -4.261; M[3,2] = -3.647; M[3,3] = -3.287 M[4,1] = -3.828; M[4,2] = -3.248; M[4,3] = -3.005 return(M) } // One break, break model: 4 sample sizes x 9 lambda = 36 rows real matrix lsmlm_cvb1() { real matrix M M = J(36, 3, .) M[1,1]=-4.630; M[1,2]=-4.064; M[1,3]=-3.787 M[2,1]=-4.704; M[2,2]=-4.132; M[2,3]=-3.843 M[3,1]=-4.769; M[3,2]=-4.199; M[3,3]=-3.906 M[4,1]=-4.820; M[4,2]=-4.253; M[4,3]=-3.963 M[5,1]=-4.857; M[5,2]=-4.293; M[5,3]=-4.008 M[6,1]=-4.891; M[6,2]=-4.324; M[6,3]=-4.042 M[7,1]=-4.899; M[7,2]=-4.338; M[7,3]=-4.058 M[8,1]=-4.910; M[8,2]=-4.348; M[8,3]=-4.071 M[9,1]=-4.915; M[9,2]=-4.351; M[9,3]=-4.073 M[10,1]=-4.523; M[10,2]=-3.974; M[10,3]=-3.697 M[11,1]=-4.536; M[11,2]=-3.982; M[11,3]=-3.702 M[12,1]=-4.563; M[12,2]=-4.010; M[12,3]=-3.726 M[13,1]=-4.609; M[13,2]=-4.060; M[13,3]=-3.777 M[14,1]=-4.647; M[14,2]=-4.098; M[14,3]=-3.819 M[15,1]=-4.667; M[15,2]=-4.125; M[15,3]=-3.846 M[16,1]=-4.672; M[16,2]=-4.130; M[16,3]=-3.856 M[17,1]=-4.681; M[17,2]=-4.144; M[17,3]=-3.870 M[18,1]=-4.685; M[18,2]=-4.145; M[18,3]=-3.873 M[19,1]=-4.502; M[19,2]=-3.959; M[19,3]=-3.679 M[20,1]=-4.489; M[20,2]=-3.944; M[20,3]=-3.664 M[21,1]=-4.496; M[21,2]=-3.950; M[21,3]=-3.671 M[22,1]=-4.534; M[22,2]=-3.991; M[22,3]=-3.714 M[23,1]=-4.570; M[23,2]=-4.031; M[23,3]=-3.756 M[24,1]=-4.587; M[24,2]=-4.052; M[24,3]=-3.777 M[25,1]=-4.598; M[25,2]=-4.065; M[25,3]=-3.794 M[26,1]=-4.612; M[26,2]=-4.081; M[26,3]=-3.809 M[27,1]=-4.602; M[27,2]=-4.074; M[27,3]=-3.799 M[28,1]=-4.466; M[28,2]=-3.928; M[28,3]=-3.651 M[29,1]=-4.455; M[29,2]=-3.919; M[29,3]=-3.640 M[30,1]=-4.469; M[30,2]=-3.933; M[30,3]=-3.651 M[31,1]=-4.501; M[31,2]=-3.968; M[31,3]=-3.687 M[32,1]=-4.523; M[32,2]=-3.993; M[32,3]=-3.713 M[33,1]=-4.546; M[33,2]=-4.018; M[33,3]=-3.740 M[34,1]=-4.562; M[34,2]=-4.033; M[34,3]=-3.756 M[35,1]=-4.576; M[35,2]=-4.049; M[35,3]=-3.773 M[36,1]=-4.593; M[36,2]=-4.065; M[36,3]=-3.790 return(M) } // Two breaks, break model: 4 sample sizes x 12 = 48 rows real matrix lsmlm_cvb2() { real matrix M M = J(48, 3, .) M[1,1]=-6.750; M[1,2]=-6.108; M[1,3]=-5.779 M[2,1]=-7.196; M[2,2]=-6.312; M[2,3]=-5.893 M[3,1]=-6.932; M[3,2]=-6.175; M[3,3]=-5.825 M[4,1]=-7.004; M[4,2]=-6.185; M[4,3]=-5.828 M[5,1]=-6.691; M[5,2]=-6.152; M[5,3]=-5.798 M[6,1]=-6.821; M[6,2]=-5.917; M[6,3]=-5.541 M[7,1]=-6.963; M[7,2]=-6.201; M[7,3]=-5.890 M[8,1]=-6.821; M[8,2]=-6.166; M[8,3]=-5.832 M[9,1]=-6.978; M[9,2]=-6.288; M[9,3]=-5.998 M[10,1]=-7.032; M[10,2]=-6.375; M[10,3]=-6.011 M[11,1]=-6.863; M[11,2]=-6.268; M[11,3]=-5.956 M[12,1]=-7.014; M[12,2]=-6.446; M[12,3]=-6.072 M[13,1]=-5.667; M[13,2]=-5.177; M[13,3]=-4.921 M[14,1]=-5.974; M[14,2]=-5.342; M[14,3]=-5.004 M[15,1]=-5.869; M[15,2]=-5.398; M[15,3]=-5.112 M[16,1]=-6.064; M[16,2]=-5.462; M[16,3]=-5.192 M[17,1]=-5.957; M[17,2]=-5.374; M[17,3]=-5.104 M[18,1]=-5.783; M[18,2]=-5.272; M[18,3]=-4.943 M[19,1]=-5.934; M[19,2]=-5.300; M[19,3]=-4.988 M[20,1]=-5.771; M[20,2]=-5.357; M[20,3]=-5.032 M[21,1]=-5.958; M[21,2]=-5.421; M[21,3]=-5.179 M[22,1]=-5.751; M[22,2]=-5.349; M[22,3]=-5.063 M[23,1]=-5.971; M[23,2]=-5.264; M[23,3]=-5.018 M[24,1]=-6.035; M[24,2]=-5.484; M[24,3]=-5.161 M[25,1]=-5.640; M[25,2]=-4.855; M[25,3]=-4.572 M[26,1]=-5.392; M[26,2]=-4.939; M[26,3]=-4.691 M[27,1]=-5.608; M[27,2]=-4.947; M[27,3]=-4.661 M[28,1]=-5.602; M[28,2]=-5.059; M[28,3]=-4.782 M[29,1]=-5.533; M[29,2]=-4.989; M[29,3]=-4.752 M[30,1]=-5.619; M[30,2]=-5.038; M[30,3]=-4.723 M[31,1]=-5.511; M[31,2]=-4.898; M[31,3]=-4.607 M[32,1]=-5.494; M[32,2]=-5.035; M[32,3]=-4.742 M[33,1]=-5.548; M[33,2]=-5.043; M[33,3]=-4.811 M[34,1]=-5.716; M[34,2]=-5.149; M[34,3]=-4.912 M[35,1]=-5.509; M[35,2]=-4.963; M[35,3]=-4.674 M[36,1]=-5.655; M[36,2]=-4.967; M[36,3]=-4.705 M[37,1]=-5.116; M[37,2]=-4.539; M[37,3]=-4.195 M[38,1]=-5.240; M[38,2]=-4.578; M[38,3]=-4.278 M[39,1]=-5.437; M[39,2]=-4.917; M[39,3]=-4.533 M[40,1]=-5.387; M[40,2]=-4.884; M[40,3]=-4.567 M[41,1]=-5.099; M[41,2]=-4.716; M[41,3]=-4.481 M[42,1]=-5.247; M[42,2]=-4.776; M[42,3]=-4.472 M[43,1]=-5.085; M[43,2]=-4.576; M[43,3]=-4.301 M[44,1]=-5.209; M[44,2]=-4.697; M[44,3]=-4.398 M[45,1]=-5.212; M[45,2]=-4.734; M[45,3]=-4.534 M[46,1]=-5.183; M[46,2]=-4.773; M[46,3]=-4.535 M[47,1]=-5.100; M[47,2]=-4.699; M[47,3]=-4.364 M[48,1]=-5.291; M[48,2]=-4.803; M[48,3]=-4.481 return(M) } // Linear interpolation over the sample size T real rowvector lsmlm_interp_T(real matrix tab, real scalar tnobs) { real rowvector nl, cv real scalar i, t1, t2, tw1, tw2 nl = lsmlm_nlookup() t1 = cols(nl); t2 = t1; tw1 = 1; tw2 = 0 for (i=1; i<=cols(nl); i++) { if (tnobs < nl[i]) { if (i == 1) { t1 = 1; t2 = 1; tw1 = 1; tw2 = 0 } else { tw1 = (tnobs - nl[i-1]) / (nl[i] - nl[i-1]) t1 = i-1 tw2 = 1 - tw1 t2 = i } break } } cv = tw1 * tab[t1,.] + tw2 * tab[t2,.] return(cv) } // One break, break model: interpolate over lambda then T real rowvector lsmlm_interp_cvb1(real scalar lambda, real scalar tnobs) { real rowvector bp, nl, cv real scalar i, i1, i2, w1, w2, lam, off real matrix M, cv1 bp = (.10, .15, .20, .25, .30, .35, .40, .45, .50) nl = lsmlm_nlookup() M = lsmlm_cvb1() lam = lambda if (lam > 0.5) lam = 1 - lam i1 = 1; i2 = 1; w1 = 1; w2 = 0 if (lam <= bp[1]) { i1 = 1; i2 = 1; w1 = 1; w2 = 0 } else if (lam >= bp[cols(bp)]) { i1 = cols(bp); i2 = i1; w1 = 1; w2 = 0 } else { for (i=2; i<=cols(bp); i++) { if (lam < bp[i]) { w2 = (lam - bp[i-1]) / (bp[i] - bp[i-1]) w1 = 1 - w2 i1 = i-1 i2 = i break } } } cv1 = J(cols(nl), 3, .) for (i=1; i<=cols(nl); i++) { off = (i-1)*9 cv1[i,.] = w1 * M[off+i1,.] + w2 * M[off+i2,.] } cv = lsmlm_interp_T(cv1, tnobs) return(cv) } // Two breaks, break model: nearest (lambda1,lambda2) cell then interpolate T real rowvector lsmlm_interp_cvb2(real scalar lam1, real scalar lam2, real scalar tnobs) { real rowvector bp, nl, cv real scalar i, i1, i2, lambda1, lambda2, temp, idx, off, b1, b2 real matrix M, cv1 bp = (.20, .30, .40, .50, .60, .70, .80) nl = lsmlm_nlookup() M = lsmlm_cvb2() lambda1 = lam1 lambda2 = lam2 if (lambda1 + lambda2 >= 1.0) { temp = lambda1 lambda1 = 1 - lambda2 lambda2 = 1 - temp } b1 = 1; b2 = 1; i1 = 1; i2 = 1 for (i=1; i<=cols(bp); i++) { if (abs(lambda1 - bp[i]) < b1) { i1 = i b1 = abs(lambda1 - bp[i]) } if (abs(lambda2 - bp[i]) < b2) { i2 = i b2 = abs(lambda2 - bp[i]) } } if (i1 == i2) { if (i1 == 1) { i2 = i1 + 1 } else if (i1 == cols(bp)) { i1 = i2 - 1 } else if ((lambda1 - bp[i1]) + (lambda2 - bp[i2]) > 0) { i2 = i1 + 1 } else { i1 = i2 - 1 } } idx = (cols(bp) - i1) * (i1 - 1) + i2 - 1 if (idx < 1) idx = 1 if (idx > 12) idx = 12 cv1 = J(cols(nl), 3, .) for (i=1; i<=cols(nl); i++) { off = (i-1)*12 cv1[i,.] = M[off+idx,.] } cv = lsmlm_interp_T(cv1, tnobs) return(cv) } // Dispatcher: pick the right table given breaks/model, write CVs to r() void lsmlm_cvlookup(real scalar breaks, real scalar model, real scalar lam1, real scalar lam2, real scalar tnobs) { real rowvector cv if (breaks == 0) { cv = lsmlm_interp_T(lsmlm_cv0(), tnobs) } else if (breaks == 1 & model == 1) { cv = lsmlm_interp_T(lsmlm_cvc1(), tnobs) } else if (breaks == 2 & model == 1) { cv = lsmlm_interp_T(lsmlm_cvc2(), tnobs) } else if (breaks == 1 & model == 2) { cv = lsmlm_interp_cvb1(lam1, tnobs) } else if (breaks == 2 & model == 2) { cv = lsmlm_interp_cvb2(lam1, lam2, tnobs) } else { cv = J(1, 3, .) } st_numscalar("r(cv1)", cv[1]) st_numscalar("r(cv5)", cv[2]) st_numscalar("r(cv10)", cv[3]) } end