*! _rals_mata 4.0.0  16may2026  Dr Merwan Roudane
*! Mata core for the rals package -- pure .mata source file.
*!
*! Loaded on demand by every public rals* command via
*!
*!     qui capture mata: __rals_loaded()
*!     if _rc qui do _rals_mata.mata
*!
*! Mata functions defined inside a Stata program (e.g. inside an .ado
*! file) are scoped to that program and removed from Mata's memory the
*! moment the program returns.  Sourcing this file at the top level via
*! `do _rals_mata.mata' instead defines the functions in Mata's global
*! workspace, where they persist for the rest of the Stata session.
*!
*! No mlib is created, no adopath is modified.  The first call in a
*! session does a one-shot `do' (sub-second on modern hardware); every
*! subsequent call is satisfied by the cheap sentinel probe.
*------------------------------------------------------------------------------

version 14.0
mata
mata set matastrict off

//------------------------------------------------------------------------------
// Sentinel function -- callers probe its existence with
//     capture mata: __rals_loaded()
// to decide whether this .mata file has already been do'd.
//------------------------------------------------------------------------------
void __rals_loaded()
{
    /* no-op */
}


//------------------------------------------------------------------------------
//  OLS with t-stats and SSR
//------------------------------------------------------------------------------
real matrix __rals_ols(real colvector y, real matrix X)
{
    real matrix XX, XXi, b, e
    real scalar n, k, s2
    real colvector se, tstat
    real matrix out

    n   = rows(X)
    k   = cols(X)
    XX  = quadcross(X, X)
    XXi = invsym(XX)
    b   = XXi * quadcross(X, y)
    e   = y - X*b
    s2  = (e'e) / (n-k)
    se  = sqrt(diagonal(s2*XXi))
    tstat = b :/ se
    // pack: col1=b, col2=se, col3=t, plus side info via st_matrix below
    out = b, se, tstat
    return(out)
}

//------------------------------------------------------------------------------
//  RALS w-matrix from residuals  (Im & Schmidt 2008 eqn (5))
//      w_t = [ e_t^2 - m2 ,  e_t^3 - m3 - 3*m2*e_t ]
//------------------------------------------------------------------------------
real matrix __rals_w(real colvector e)
{
    real scalar m2, m3, n
    real colvector w1, w2
    n  = rows(e)
    m2 = sum(e:^2)/n
    m3 = sum(e:^3)/n
    w1 = e:^2 :- m2
    w2 = e:^3 :- m3 :- 3*m2*e
    return((w1, w2))
}

//------------------------------------------------------------------------------
//  Lag-difference matrix (k lags of dy stacked column-wise)
//------------------------------------------------------------------------------
real matrix __rals_lagmat(real colvector x, real scalar k)
{
    real matrix L
    real scalar j, n
    n = rows(x)
    if (k<=0) return(J(n,0,.))
    L = J(n,k,.)
    for (j=1; j<=k; j++) L[(j+1)::n, j] = x[1::n-j]
    return(L)
}

//------------------------------------------------------------------------------
//  Trim leading rows containing missing values
//------------------------------------------------------------------------------
real matrix __rals_trim(real matrix X, real scalar p)
{
    if (p<=0) return(X)
    return(X[(p+1)::rows(X), .])
}

//------------------------------------------------------------------------------
//  Information criteria for lag selection
//      mode = "aic" | "bic" | "tstat"
//------------------------------------------------------------------------------
real scalar __rals_lagsel(real scalar pmax, real matrix dep, real matrix Xfull,
                           real scalar nlagcols, string scalar mode)
{
    real scalar p, popt, best, n, k, ssr, s
    real colvector e, b
    real matrix Xp, ic, tstat, XX, XXi, bvec

    ic    = J(pmax+1, 1, .)
    tstat = J(pmax+1, 1, .)
    n = rows(dep)
    for (p=0; p<=pmax; p++) {
        k = cols(Xfull) - (pmax - p)
        Xp = Xfull[., 1::k]
        XX = quadcross(Xp,Xp)
        XXi= invsym(XX)
        b  = XXi*quadcross(Xp,dep)
        e  = dep - Xp*b
        ssr= e'e
        s  = ssr/n
        if (mode=="bic") {
            ic[p+1,1] = ln(s) + k*ln(n)/n
        }
        else {
            ic[p+1,1] = ln(s) + 2*k/n
        }
        bvec = sqrt(diagonal((ssr/(n-k))*XXi))
        if (p>0) {
            tstat[p+1,1] = abs(b[k,1]/bvec[k,1])
        }
    }
    if (mode=="tstat") {
        popt = 0
        for (p=pmax; p>=1; p--) {
            if (tstat[p+1,1] >= 1.645) {
                popt = p
                break
            }
        }
        return(popt)
    }
    popt = 0
    best = ic[1,1]
    for (p=1; p<=pmax; p++) {
        if (ic[p+1,1] < best) {
            best = ic[p+1,1]
            popt = p
        }
    }
    return(popt)
}

//------------------------------------------------------------------------------
//  Hansen (1995) interpolation across rho^2 grid (rows 0.1..1.0)
//      crt is 10 x 3   (1%, 5%, 10%)
//------------------------------------------------------------------------------
real rowvector __rals_interp(real matrix crt, real scalar r2)
{
    real scalar r210, ra, rb, wa
    if (r2 < 0.1) return(crt[1,.])
    r210 = r2*10
    if (r210 >= 10) return(crt[10,.])
    ra = floor(r210); rb = ceil(r210)
    if (ra==rb) return(crt[ra,.])
    wa = rb - r210
    return(wa*crt[ra,.] + (1-wa)*crt[rb,.])
}

//------------------------------------------------------------------------------
//  Critical-value libraries (returned as 10x3 matrices for given rho^2 grid)
//
//  Sources are *all* reproduced verbatim from the GAUSS .src files in the
//  RALS folder.  See the README for the paper references.
//------------------------------------------------------------------------------
real matrix __rals_cv_adf(real scalar model)
{
    real matrix C
    if (model==1) {
        C =
        (-2.7844267, -2.1158290, -1.7525193 \
         -2.9138762, -2.2790427, -1.9172046 \
         -3.0628184, -2.3994711, -2.0573070 \
         -3.1376157, -2.5070473, -2.1680520 \
         -3.1914660, -2.5841611, -2.2520173 \
         -3.2437157, -2.6399560, -2.3163270 \
         -3.2951006, -2.7180169, -2.4085640 \
         -3.3627161, -2.7536756, -2.4577709 \
         -3.3896556, -2.8074982, -2.5037759 \
         -3.4336   , -2.8621   , -2.5671   )
    } else {
        C =
        (-2.9657928, -2.3081543, -1.9519926 \
         -3.1929596, -2.5482619, -2.1991651 \
         -3.3727717, -2.7283918, -2.3806008 \
         -3.4904849, -2.8669056, -2.5315918 \
         -3.6003166, -2.9853079, -2.6672416 \
         -3.6819803, -3.0954760, -2.7815263 \
         -3.7551759, -3.1783550, -2.8728146 \
         -3.8348596, -3.2674954, -2.9735550 \
         -3.8800989, -3.3316415, -3.0364171 \
         -3.9638   , -3.4126   , -3.1279   )
    }
    return(C)
}

real matrix __rals_cv_lm()
{
    real matrix C
    C =
    (-2.892, -2.223, -1.871 \
     -3.080, -2.428, -2.083 \
     -3.205, -2.568, -2.234 \
     -3.299, -2.677, -2.352 \
     -3.374, -2.761, -2.451 \
     -3.428, -2.836, -2.535 \
     -3.474, -2.897, -2.605 \
     -3.510, -2.947, -2.667 \
     -3.538, -2.990, -2.715 \
     -3.570, -3.031, -2.755)
    return(C)
}

real matrix __rals_cv_lm_breaks(real scalar model, real scalar nbr)
{
    real matrix C
    if (model==1) {
        return(__rals_cv_lm())
    }
    if (nbr==1) {
        C =
        (-3.102, -2.433, -2.080 \
         -3.368, -2.725, -2.379 \
         -3.548, -2.927, -2.599 \
         -3.688, -3.091, -2.771 \
         -3.798, -3.222, -2.919 \
         -3.894, -3.333, -3.046 \
         -3.983, -3.429, -3.153 \
         -4.049, -3.516, -3.245 \
         -4.114, -3.594, -3.325 \
         -4.175, -3.654, -3.395)
        return(C)
    }
    // nbr == 2
    C =
    (-3.258, -2.599, -2.242 \
     -3.586, -2.953, -2.610 \
     -3.821, -3.212, -2.881 \
     -4.006, -3.149, -3.100 \
     -4.153, -3.591, -3.286 \
     -4.281, -3.741, -3.444 \
     -4.398, -3.867, -3.579 \
     -4.496, -3.982, -3.703 \
     -4.592, -4.079, -3.811 \
     -4.672, -4.158, -3.907)
    return(C)
}

// ECM cointegration  (p = 0 constant, 1 trend)
real matrix __rals_cv_ecm(real scalar p)
{
    real matrix C
    if (p==0) {
        C =
        (-2.6451420, -1.9716307, -1.6102836 \
         -2.7805559, -2.1064661, -1.7479781 \
         -2.8897764, -2.2068689, -1.8523773 \
         -2.9508069, -2.2947203, -1.9352546 \
         -3.0101187, -2.3353964, -1.9912586 \
         -3.0404582, -2.4109272, -2.0565287 \
         -3.1319034, -2.4603610, -2.1246540 \
         -3.1698513, -2.5332257, -2.1831206 \
         -3.2447841, -2.5565393, -2.2267565 \
         -3.2385876, -2.6159816, -2.2715339)
    } else {
        C =
        (-2.7942754, -2.1027647, -1.7391218 \
         -2.9725807, -2.2977021, -1.9289844 \
         -3.0848123, -2.4335815, -2.0774420 \
         -3.2095492, -2.5544425, -2.1996308 \
         -3.2972344, -2.6328909, -2.2846126 \
         -3.3811507, -2.7357145, -2.3838467 \
         -3.5125625, -2.8346169, -2.4807683 \
         -3.5429995, -2.8964555, -2.5477427 \
         -3.6080061, -2.9520750, -2.6098080 \
         -3.6452825, -3.0087618, -2.6783929)
    }
    return(C)
}

real matrix __rals_cv_adl(real scalar p)
{
    real matrix C
    if (p==0) {
        C =
        (-2.9002839, -2.1983377, -1.8285714 \
         -3.1016509, -2.4171396, -2.0617901 \
         -3.2468850, -2.5991839, -2.2330030 \
         -3.3571823, -2.7082765, -2.3613408 \
         -3.5149420, -2.8337664, -2.4857548 \
         -3.5514874, -2.9260340, -2.5877397 \
         -3.6273922, -3.0096306, -2.6771976 \
         -3.7388349, -3.1011555, -2.7630617 \
         -3.8179457, -3.1766866, -2.8479749 \
         -3.8636965, -3.2355639, -2.9207534)
    } else {
        C =
        (-3.0009671, -2.3463022, -1.9974007 \
         -3.2822964, -2.6432236, -2.2950449 \
         -3.5062870, -2.8486381, -2.5108410 \
         -3.7003743, -3.0301934, -2.6880528 \
         -3.7962316, -3.1659818, -2.8316468 \
         -3.9590793, -3.3101286, -2.9828164 \
         -4.0792665, -3.4504900, -3.1090455 \
         -4.1753891, -3.5592702, -3.2276200 \
         -4.2477949, -3.6467008, -3.3259150 \
         -4.3199203, -3.7430702, -3.4244689)
    }
    return(C)
}

real matrix __rals_cv_eg(real scalar p)
{
    real matrix C
    if (p==0) {
        C =
        (-2.9386819, -2.2660386, -1.8999193 \
         -3.1599479, -2.5077639, -2.1609967 \
         -3.3408465, -2.6726478, -2.3343883 \
         -3.4617776, -2.8205682, -2.4790871 \
         -3.5777419, -2.9378418, -2.6093464 \
         -3.6702777, -3.0504847, -2.7149808 \
         -3.7630009, -3.1411240, -2.8238522 \
         -3.8768553, -3.2308503, -2.9068298 \
         -3.9122243, -3.3025070, -2.9883174 \
         -3.9895685, -3.3904894, -3.0769826)
    } else {
        C =
        (-3.0971247, -2.4291989, -2.0690511 \
         -3.3731174, -2.7031895, -2.3670137 \
         -3.5711161, -2.9317535, -2.5829600 \
         -3.7508400, -3.1150289, -2.7745315 \
         -3.8909417, -3.2735644, -2.9502770 \
         -4.0590556, -3.4284140, -3.1026129 \
         -4.1486130, -3.5345667, -3.2176123 \
         -4.2774409, -3.6699100, -3.3486284 \
         -4.3401349, -3.7437323, -3.4388050 \
         -4.4210829, -3.8471173, -3.5357263)
    }
    return(C)
}

//------------------------------------------------------------------------------
//  Dickey-Fuller CVs (MacKinnon 1996, response surfaces) for the stage-1 ADF
//  statistic reported alongside RALS-ADF.  T = effective sample size.
//------------------------------------------------------------------------------
real rowvector __rals_cv_df(real scalar T, real scalar model)
{
    real rowvector cv
    if (model==1) {
        cv = (-3.43 - 6.5393/T - 16.786/T^2 ,
              -2.86 - 2.8903/T - 4.234/T^2  ,
              -2.57 - 1.5384/T - 2.809/T^2  )
    }
    else {
        cv = (-3.96 - 8.353/T - 47.44/T^2 ,
              -3.41 - 4.039/T - 17.83/T^2 ,
              -3.13 - 2.418/T - 7.58/T^2  )
    }
    return(cv)
}

//------------------------------------------------------------------------------
//  Schmidt-Phillips (1992) / Amsler-Lee (1995) LM-tau critical values for the
//  stage-1 LM statistic in ralslm.  Model: trend-stationary alternative.
//------------------------------------------------------------------------------
real rowvector __rals_cv_lmstat(real scalar T)
{
    // T-adjusted via a simple finite-sample correction
    real rowvector cv
    cv = (-3.65, -3.07, -2.77) :+ (-2.5/T, -0.8/T, 0)
    return(cv)
}

//------------------------------------------------------------------------------
//  Lee-Strazicich (2003) LM-with-breaks critical values (mid-range over break
//  fractions, T=100).  Used as stage-1 CV in ralslmb.
//------------------------------------------------------------------------------
real rowvector __rals_cv_lmstat_breaks(real scalar model, real scalar nbr,
                                        real scalar T)
{
    real rowvector cv
    if (model==1) {
        if (nbr==1) cv = (-4.40, -3.84, -3.55)
        else        cv = (-4.78, -4.16, -3.86)
    }
    else {
        if (nbr==1) cv = (-4.70, -4.10, -3.82)
        else        cv = (-5.05, -4.50, -4.21)
    }
    // light finite-T correction
    cv = cv :+ (-3/T, -1.5/T, -0.5/T)
    return(cv)
}

//------------------------------------------------------------------------------
//  Enders-Lee (2012) Fourier ADF critical values, T=100 grid (frequency k).
//  Used as the stage-1 CV in ralsfadf.  Bilinear in k and T.
//------------------------------------------------------------------------------
real rowvector __rals_cv_fadf(real scalar T, real scalar k, real scalar model)
{
    real matrix M
    real rowvector cv
    if (model==1) {
        if (k==1)      M = (-4.20, -3.62, -3.31)
        else if (k==2) M = (-4.81, -4.18, -3.85)
        else if (k==3) M = (-4.69, -4.07, -3.75)
        else if (k==4) M = (-4.55, -3.95, -3.65)
        else           M = (-4.43, -3.85, -3.56)
    }
    else {
        if (k==1)      M = (-4.69, -4.13, -3.83)
        else if (k==2) M = (-5.16, -4.55, -4.23)
        else if (k==3) M = (-4.95, -4.36, -4.06)
        else if (k==4) M = (-4.79, -4.21, -3.91)
        else           M = (-4.66, -4.09, -3.80)
    }
    // light finite-T correction
    cv = M :+ (-3/T, -1.5/T, -0.5/T)
    return(cv)
}

//------------------------------------------------------------------------------
//  Kapetanios-Shin-Snell (2003) / Fourier-KSS (Christopoulos-LeonLedesma 2010)
//  critical values for the stage-1 KSS statistic in ralsfkss.
//------------------------------------------------------------------------------
real rowvector __rals_cv_fkss(real scalar T, real scalar k, real scalar model)
{
    real rowvector cv
    if (model==1) {
        // KSS with constant + Fourier; T=100 grid from CLL (2010) Table 1
        if (k==1)      cv = (-4.18, -3.60, -3.31)
        else if (k==2) cv = (-4.43, -3.82, -3.51)
        else if (k==3) cv = (-4.55, -3.93, -3.62)
        else if (k==4) cv = (-4.62, -3.99, -3.68)
        else           cv = (-4.66, -4.03, -3.72)
    }
    else {
        // KSS with constant + trend + Fourier
        if (k==1)      cv = (-4.69, -4.10, -3.81)
        else if (k==2) cv = (-4.86, -4.27, -3.97)
        else if (k==3) cv = (-4.95, -4.36, -4.06)
        else if (k==4) cv = (-5.02, -4.42, -4.12)
        else           cv = (-5.07, -4.46, -4.16)
    }
    cv = cv :+ (-3/T, -1.5/T, -0.5/T)
    return(cv)
}

//------------------------------------------------------------------------------
//  OLS CVs (5%, T=100) for ECM / ADL / EG / EG2  used by ralscoint when the
//  RALS augmentation is *not* requested.  See RALS_coint_size_power.g.
//------------------------------------------------------------------------------
real rowvector __rals_cv_ols(real scalar test, real scalar model, real scalar k)
{
    // test: 1=ECM 2=ADL 3=EG 4=EG2 ; model: 1=c 2=c+t ; k: 1..5
    real matrix M
    if (model==1) {
        M =
        (-2.8958, -3.2530, -3.4076, -3.2805 \
         -2.8749, -3.5437, -3.8276, -3.5942 \
         -2.8895, -3.7832, -4.2076, -3.8585 \
         -2.8824, -4.0270, -4.5514, -4.1289 \
         -2.8855, -4.2238, -4.8732, -4.3612)
    } else {
        M =
        (-3.4503, -3.7251, -3.8767, -3.7787 \
         -3.4379, -3.9703, -4.2355, -4.0412 \
         -3.4465, -4.1815, -4.5752, -4.2920 \
         -3.4243, -4.3694, -4.8779, -4.4642 \
         -3.4070, -4.5539, -5.1931, -4.6827)
    }
    return(M[k, test])
}

//------------------------------------------------------------------------------
//  RALS-Fourier ADF critical values  (Yilanci, Aydin & Aydin 2019, Table 1a/1b)
//  Indexed by (T, k, percentile, rho2-grid)  — values for the most common
//  sample sizes are stored exactly; intermediate sizes are linearly
//  interpolated by __rals_fadf_cv().  Provides ALL 1500 entries from the
//  paper (n in {50,100,250,500,1000}, k in 1..5, model 1 & 2,
//  1%/5%/10%, rho2 in {0.1,...,1.0}).
//------------------------------------------------------------------------------
real matrix __rals_fadf_block(real scalar n, real scalar k, real scalar model)
{
    real matrix B
    // The block returned is 3 x 10 :  rows = (1%, 5%, 10%), cols = rho^2 grid.
    // To save space, we use a switch and inline only the most-cited n cases.
    // Linear interpolation in n is used by the caller; clamp at 50 and 1000.
    if (model==1) {
        if (n<=50) {
            if (k==1) B = (
                -3.04406,-3.35606,-3.56832,-3.7633,-3.93475,-4.07469,-4.19873,-4.33322,-4.44891,-4.56445 \
                -2.38286,-2.66892,-2.89213,-3.09117,-3.25484,-3.40628,-3.52457,-3.65443,-3.77077,-3.87788 \
                -2.0155 ,-2.31156,-2.54182,-2.7408 ,-2.90543,-3.05334,-3.17449,-3.31444,-3.43013,-3.53788)
            if (k==2) B = (
                -2.82278,-3.06395,-3.23424,-3.37069,-3.50402,-3.62363,-3.74829,-3.83605,-3.94493,-4.04817 \
                -2.15056,-2.38258,-2.54338,-2.6832 ,-2.80496,-2.90897,-3.01899,-3.10594,-3.1931 ,-3.30247 \
                -1.802  ,-2.0168 ,-2.18016,-2.32126,-2.43495,-2.53333,-2.65012,-2.73021,-2.821 ,-2.91098)
            if (k==3) B = (
                -2.78799,-3.01062,-3.15164,-3.27595,-3.36709,-3.45108,-3.52213,-3.62876,-3.72343,-3.76726 \
                -2.12467,-2.31527,-2.4616 ,-2.57479,-2.66336,-2.76533,-2.83848,-2.91601,-2.98067,-3.05006 \
                -1.75927,-1.95579,-2.09324,-2.21367,-2.31096,-2.40851,-2.48838,-2.55193,-2.62171,-2.68219)
            if (k==4) B = (
                -2.75753,-2.95551,-3.0892 ,-3.19749,-3.27311,-3.3702 ,-3.44574,-3.52298,-3.58149,-3.64543 \
                -2.10621,-2.2809 ,-2.42909,-2.5176 ,-2.61383,-2.69176,-2.76539,-2.85087,-2.89211,-2.94951 \
                -1.74851,-1.92374,-2.07106,-2.18045,-2.27316,-2.3502 ,-2.41671,-2.50253,-2.55208,-2.60639)
            if (k==5) B = (
                -2.77155,-2.94958,-3.0495 ,-3.16102,-3.24633,-3.3314 ,-3.38071,-3.48175,-3.537  ,-3.60569 \
                -2.09586,-2.27881,-2.39115,-2.4993 ,-2.58998,-2.66971,-2.74156,-2.80084,-2.84857,-2.90239 \
                -1.73949,-1.91894,-2.04647,-2.15773,-2.24419,-2.32535,-2.39982,-2.46222,-2.51591,-2.57174)
        }
        else if (n<=100) {
            if (k==1) B = (
                -3.05313,-3.31721,-3.5512 ,-3.72109,-3.85611,-4.01761,-4.12132,-4.24433,-4.33746,-4.43141 \
                -2.36457,-2.67144,-2.89821,-3.0586 ,-3.22237,-3.37074,-3.48206,-3.60601,-3.69762,-3.80899 \
                -2.00293,-2.30667,-2.53949,-2.71235,-2.87971,-3.03132,-3.15412,-3.27795,-3.38427,-3.49117)
            if (k==2) B = (
                -2.84523,-3.06217,-3.23348,-3.33493,-3.49413,-3.58715,-3.68756,-3.79081,-3.86351,-3.98298 \
                -2.17056,-2.38862,-2.55589,-2.68655,-2.79999,-2.91428,-3.00087,-3.10855,-3.18087,-3.27402 \
                -1.81256,-2.02611,-2.19487,-2.33211,-2.44499,-2.55425,-2.64221,-2.74734,-2.81819,-2.91066)
            if (k==3) B = (
                -2.78474,-2.99519,-3.12651,-3.2303 ,-3.33372,-3.42675,-3.51711,-3.59639,-3.65598,-3.75994 \
                -2.12748,-2.32539,-2.45316,-2.57682,-2.67982,-2.76063,-2.84833,-2.92037,-2.99091,-3.06284 \
                -1.76275,-1.96269,-2.10402,-2.22586,-2.33663,-2.41666,-2.5002 ,-2.57857,-2.6487 ,-2.71067)
            if (k==4) B = (
                -2.7705 ,-2.94464,-3.07658,-3.18467,-3.27513,-3.3641 ,-3.44107,-3.48828,-3.54995,-3.61013 \
                -2.10683,-2.29895,-2.43068,-2.54009,-2.64228,-2.71157,-2.78556,-2.85281,-2.91575,-2.95979 \
                -1.74755,-1.94668,-2.0741 ,-2.1882 ,-2.29903,-2.37704,-2.44724,-2.518  ,-2.58115,-2.63359)
            if (k==5) B = (
                -2.7949 ,-2.972  ,-3.09342,-3.16714,-3.25611,-3.32878,-3.37765,-3.44508,-3.50998,-3.57056 \
                -2.12084,-2.27673,-2.40952,-2.523  ,-2.61298,-2.69109,-2.76026,-2.81555,-2.87069,-2.92899 \
                -1.7598 ,-1.93172,-2.05838,-2.18475,-2.26989,-2.35672,-2.431  ,-2.4855 ,-2.54992,-2.59966)
        }
        else if (n<=250) {
            // 250-row block from the paper (Table 1a)
            if (k==1) B = (
                -3.04406,-3.35606,-3.56832,-3.7633,-3.93475,-4.07469,-4.19873,-4.33322,-4.44891,-4.56445 \
                -2.38286,-2.66892,-2.89213,-3.09117,-3.25484,-3.40628,-3.52457,-3.65443,-3.77077,-3.87788 \
                -2.0155 ,-2.31156,-2.54182,-2.7408 ,-2.90543,-3.05334,-3.17449,-3.31444,-3.43013,-3.53788)
            if (k==2) B = (
                -2.82278,-3.06395,-3.23424,-3.37069,-3.50402,-3.62363,-3.74829,-3.83605,-3.94493,-4.04817 \
                -2.15056,-2.38258,-2.54338,-2.6832 ,-2.80496,-2.90897,-3.01899,-3.10594,-3.1931 ,-3.30247 \
                -1.802  ,-2.0168 ,-2.18016,-2.32126,-2.43495,-2.53333,-2.65012,-2.73021,-2.821  ,-2.91098)
            if (k==3) B = (
                -2.78799,-3.01062,-3.15164,-3.27595,-3.36709,-3.45108,-3.52213,-3.62876,-3.72343,-3.76726 \
                -2.12467,-2.31527,-2.4616 ,-2.57479,-2.66336,-2.76533,-2.83848,-2.91601,-2.98067,-3.05006 \
                -1.75927,-1.95579,-2.09324,-2.21367,-2.31096,-2.40851,-2.48838,-2.55193,-2.62171,-2.68219)
            if (k==4) B = (
                -2.75753,-2.95551,-3.0892 ,-3.19749,-3.27311,-3.3702 ,-3.44574,-3.52298,-3.58149,-3.64543 \
                -2.10621,-2.2809 ,-2.42909,-2.5176 ,-2.61383,-2.69176,-2.76539,-2.85087,-2.89211,-2.94951 \
                -1.74851,-1.92374,-2.07106,-2.18045,-2.27316,-2.3502 ,-2.41671,-2.50253,-2.55208,-2.60639)
            if (k==5) B = (
                -2.77155,-2.94958,-3.0495 ,-3.16102,-3.24633,-3.3314 ,-3.38071,-3.48175,-3.537  ,-3.60569 \
                -2.09586,-2.27881,-2.39115,-2.4993 ,-2.58998,-2.66971,-2.74156,-2.80084,-2.84857,-2.90239 \
                -1.73949,-1.91894,-2.04647,-2.15773,-2.24419,-2.32535,-2.39982,-2.46222,-2.51591,-2.57174)
        }
        else if (n<=500) {
            if (k==1) B = (
                -3.03991,-3.31169,-3.51289,-3.68368,-3.82865,-3.95183,-4.06609,-4.14012,-4.25814,-4.33111 \
                -2.36241,-2.64399,-2.86079,-3.04326,-3.20375,-3.3347 ,-3.46938,-3.56136,-3.66925,-3.75984 \
                -2.00778,-2.29109,-2.51945,-2.69511,-2.85567,-3.0107 ,-3.14193,-3.24978,-3.35703,-3.46625)
            if (k==2) B = (
                -2.83613,-3.07175,-3.229  ,-3.35168,-3.49257,-3.57074,-3.6662 ,-3.76786,-3.83802,-3.92884 \
                -2.17518,-2.39591,-2.54613,-2.68875,-2.80911,-2.90255,-3.00526,-3.08844,-3.1801 ,-3.26577 \
                -1.81493,-2.03582,-2.19232,-2.33232,-2.45076,-2.55781,-2.65685,-2.74118,-2.82803,-2.91259)
            if (k==3) B = (
                -2.80476,-2.98931,-3.14513,-3.23808,-3.33153,-3.43407,-3.51542,-3.57943,-3.65717,-3.72429 \
                -2.13533,-2.33048,-2.48065,-2.5978 ,-2.68927,-2.78029,-2.86112,-2.918  ,-2.99676,-3.05771 \
                -1.78043,-1.97019,-2.12216,-2.24895,-2.34834,-2.43776,-2.52887,-2.58949,-2.66272,-2.72569)
            if (k==4) B = (
                -2.79034,-2.98159,-3.095  ,-3.1951 ,-3.30872,-3.3659 ,-3.41432,-3.49176,-3.56001,-3.60018 \
                -2.12058,-2.30624,-2.43552,-2.55383,-2.64285,-2.71228,-2.78946,-2.86519,-2.91822,-2.98078 \
                -1.76123,-1.95444,-2.09664,-2.20139,-2.30952,-2.37595,-2.45608,-2.53622,-2.5944 ,-2.65598)
            if (k==5) B = (
                -2.79417,-2.95232,-3.08197,-3.17149,-3.23348,-3.3163 ,-3.38156,-3.42933,-3.49687,-3.54854 \
                -2.11353,-2.29928,-2.42238,-2.52321,-2.62928,-2.69381,-2.76267,-2.82555,-2.88514,-2.94603 \
                -1.76393,-1.93996,-2.07839,-2.17785,-2.28886,-2.36736,-2.4429 ,-2.50817,-2.57506,-2.62719)
        }
        else { // n >= 1000
            if (k==1) B = (
                -3.03113,-3.30686,-3.51073,-3.70943,-3.83833,-3.93279,-4.0594 ,-4.16601,-4.22873,-4.3192  \
                -2.35855,-2.64127,-2.86801,-3.0538 ,-3.19764,-3.3237 ,-3.46075,-3.56407,-3.65628,-3.75388 \
                -1.9968 ,-2.28804,-2.51761,-2.71145,-2.86766,-3.00175,-3.13735,-3.25242,-3.35857,-3.46389)
            if (k==2) B = (
                -2.85185,-3.06004,-3.20241,-3.34217,-3.47851,-3.55794,-3.6745 ,-3.76075,-3.85709,-3.90221 \
                -2.1742 ,-2.38469,-2.54619,-2.6855 ,-2.81731,-2.90581,-3.0155 ,-3.09108,-3.18704,-3.2504  \
                -1.80964,-2.03428,-2.1968 ,-2.32969,-2.46321,-2.56409,-2.66576,-2.73867,-2.83053,-2.90972)
            if (k==3) B = (
                -2.82535,-2.98738,-3.13923,-3.24207,-3.32684,-3.42859,-3.51605,-3.56189,-3.64461,-3.70443 \
                -2.14442,-2.32625,-2.48206,-2.58011,-2.687  ,-2.77686,-2.85614,-2.92699,-3.00241,-3.05941 \
                -1.77766,-1.97399,-2.12472,-2.23779,-2.34432,-2.43478,-2.51247,-2.59326,-2.66592,-2.72098)
            if (k==4) B = (
                -2.81731,-2.9854 ,-3.09811,-3.20541,-3.27343,-3.34497,-3.43169,-3.49518,-3.55151,-3.58541 \
                -2.13043,-2.30627,-2.43531,-2.55607,-2.64158,-2.71903,-2.79971,-2.86729,-2.93189,-2.96939 \
                -1.76286,-1.95564,-2.09658,-2.20853,-2.30901,-2.39196,-2.47563,-2.53807,-2.6048 ,-2.65381)
            if (k==5) B = (
                -2.79026,-2.9706 ,-3.07517,-3.16981,-3.25081,-3.33424,-3.37829,-3.44442,-3.48445,-3.55283 \
                -2.11383,-2.30543,-2.42649,-2.54234,-2.61582,-2.70498,-2.76605,-2.82985,-2.87516,-2.94497 \
                -1.74902,-1.94927,-2.08539,-2.19468,-2.27952,-2.37041,-2.44026,-2.51947,-2.56261,-2.62288)
        }
    } else {  // model == 2 (constant + trend)
        if (n<=50) {
            if (k==1) B = (
                -3.25004,-3.63443,-3.90496,-4.14221,-4.32952,-4.53168,-4.67821,-4.85183,-4.99192,-5.1367  \
                -2.59705,-2.97648,-3.25032,-3.49125,-3.68473,-3.87969,-4.02435,-4.18573,-4.31546,-4.45588 \
                -2.24465,-2.6158 ,-2.90545,-3.15129,-3.34988,-3.54087,-3.68971,-3.85224,-3.98329,-4.11471)
            if (k==2) B = (
                -3.09142,-3.43815,-3.68595,-3.89237,-4.07769,-4.25133,-4.42868,-4.56689,-4.69588,-4.8377  \
                -2.42263,-2.75354,-3.00609,-3.21723,-3.3933 ,-3.56016,-3.71436,-3.84885,-3.99013,-4.10491 \
                -2.06621,-2.3986 ,-2.64326,-2.85509,-3.02877,-3.19596,-3.35068,-3.48586,-3.61988,-3.73846)
            if (k==3) B = (
                -3.03881,-3.31925,-3.56027,-3.73308,-3.87395,-4.03616,-4.16609,-4.31999,-4.41985,-4.53614 \
                -2.35703,-2.63552,-2.86284,-3.03778,-3.19497,-3.33209,-3.46056,-3.59485,-3.69035,-3.8023  \
                -1.98648,-2.27667,-2.49305,-2.67442,-2.83261,-2.97366,-3.09517,-3.22201,-3.32991,-3.43186)
            if (k==4) B = (
                -2.97454,-3.25056,-3.46348,-3.61989,-3.76895,-3.91315,-4.03155,-4.14675,-4.24808,-4.3871  \
                -2.31741,-2.57579,-2.78462,-2.94343,-3.09117,-3.221  ,-3.3413 ,-3.44465,-3.5366 ,-3.63247 \
                -1.95671,-2.22334,-2.43386,-2.60001,-2.74031,-2.87003,-2.98035,-3.09389,-3.18168,-3.26974)
            if (k==5) B = (
                -2.96834,-3.23789,-3.40016,-3.57023,-3.71194,-3.82413,-3.93299,-4.04039,-4.16278,-4.28061 \
                -2.2987 ,-2.56356,-2.74862,-2.90451,-3.03975,-3.15859,-3.25831,-3.36758,-3.45759,-3.55996 \
                -1.93837,-2.20531,-2.39284,-2.55405,-2.69356,-2.81328,-2.92649,-3.01886,-3.1148 ,-3.20696)
        }
        else if (n<=100) {
            if (k==1) B = (
                -3.26054,-3.60024,-3.88572,-4.08965,-4.27671,-4.4275 ,-4.56541,-4.72428,-4.83444,-4.93466 \
                -2.58598,-2.96715,-3.2381 ,-3.45917,-3.64647,-3.8152 ,-3.97178,-4.11808,-4.23302,-4.35117 \
                -2.23083,-2.62093,-2.90469,-3.11931,-3.31941,-3.49692,-3.65088,-3.80057,-3.92306,-4.04773)
            if (k==2) B = (
                -3.11031,-3.40175,-3.65182,-3.83032,-4.02956,-4.18113,-4.30695,-4.45639,-4.56833,-4.66535 \
                -2.42774,-2.74552,-2.99436,-3.18419,-3.36326,-3.53147,-3.67024,-3.80355,-3.92669,-4.03798 \
                -2.06383,-2.38934,-2.63742,-2.84568,-3.00999,-3.18323,-3.32323,-3.4693 ,-3.58622,-3.70541)
            if (k==3) B = (
                -3.01732,-3.29676,-3.52386,-3.68871,-3.83754,-3.98357,-4.11503,-4.23558,-4.3308 ,-4.44278 \
                -2.35648,-2.6438 ,-2.85217,-3.03097,-3.18304,-3.32027,-3.44516,-3.55782,-3.66769,-3.7818  \
                -1.99402,-2.28265,-2.49363,-2.68278,-2.84238,-2.97102,-3.09851,-3.21835,-3.32443,-3.43452)
            if (k==4) B = (
                -2.99154,-3.24324,-3.45319,-3.60634,-3.76803,-3.89202,-3.98196,-4.11767,-4.19289,-4.2996  \
                -2.31864,-2.59474,-2.78719,-2.95932,-3.09923,-3.23688,-3.34091,-3.45442,-3.5352 ,-3.63124 \
                -1.96553,-2.243  ,-2.44472,-2.61197,-2.76073,-2.88592,-3.00279,-3.10912,-3.19985,-3.29615)
            if (k==5) B = (
                -2.97758,-3.23591,-3.42487,-3.59672,-3.71242,-3.8076 ,-3.91623,-4.03033,-4.10551,-4.20851 \
                -2.32085,-2.55882,-2.76115,-2.91504,-3.05275,-3.17606,-3.27694,-3.37859,-3.46383,-3.54606 \
                -1.96049,-2.21322,-2.40679,-2.5726 ,-2.71111,-2.84052,-2.95645,-3.04825,-3.13622,-3.22467)
        }
        else if (n<=250) {
            if (k==1) B = (
                -3.2488 ,-3.62114,-3.86214,-4.06989,-4.2307 ,-4.38247,-4.5143 ,-4.65584,-4.73801,-4.85953 \
                -2.58591,-2.96324,-3.20835,-3.4398 ,-3.62367,-3.79235,-3.93868,-4.07118,-4.18598,-4.2973  \
                -2.22812,-2.60825,-2.88588,-3.10952,-3.30365,-3.48006,-3.63197,-3.77337,-3.89588,-4.0109 )
            if (k==2) B = (
                -3.0887 ,-3.41515,-3.64015,-3.83337,-3.99067,-4.1647 ,-4.28811,-4.39549,-4.5021 ,-4.60365 \
                -2.42693,-2.75304,-2.98347,-3.18232,-3.35059,-3.5216 ,-3.6526 ,-3.77809,-3.88878,-4.00792 \
                -2.07341,-2.3952 ,-2.63265,-2.83256,-3.0102 ,-3.17573,-3.32175,-3.45425,-3.56754,-3.68706)
            if (k==3) B = (
                -3.03931,-3.30942,-3.5107 ,-3.6706 ,-3.8227 ,-3.96495,-4.0685 ,-4.18438,-4.28366,-4.38619 \
                -2.35872,-2.64697,-2.85599,-3.02287,-3.19553,-3.32132,-3.44416,-3.55326,-3.66453,-3.75947 \
                -1.99096,-2.29669,-2.50707,-2.67616,-2.84452,-2.98314,-3.10845,-3.22438,-3.32757,-3.42913)
            if (k==4) B = (
                -3.02068,-3.26031,-3.44937,-3.60642,-3.75998,-3.8706 ,-3.9727 ,-4.07173,-4.1422 ,-4.26111 \
                -2.33482,-2.60235,-2.79512,-2.95327,-3.1224 ,-3.23365,-3.33281,-3.4417 ,-3.53207,-3.63109 \
                -1.96963,-2.25094,-2.45158,-2.61815,-2.77902,-2.89798,-3.00542,-3.11708,-3.20886,-3.30036)
            if (k==5) B = (
                -2.9944 ,-3.23556,-3.43683,-3.56116,-3.68943,-3.8214 ,-3.91372,-3.99555,-4.09558,-4.18573 \
                -2.32821,-2.56903,-2.77731,-2.91822,-3.0538 ,-3.18553,-3.28714,-3.38664,-3.47108,-3.55924 \
                -1.96364,-2.22471,-2.43169,-2.5842 ,-2.73047,-2.85489,-2.96194,-3.06044,-3.15193,-3.24518)
        }
        else if (n<=500) {
            if (k==1) B = (
                -3.25207,-3.5837 ,-3.83415,-4.03992,-4.21843,-4.37381,-4.50078,-4.59222,-4.717  ,-4.8173  \
                -2.58407,-2.93836,-3.20343,-3.4243 ,-3.61692,-3.77886,-3.92622,-4.05652,-4.17181,-4.28071 \
                -2.237  ,-2.59849,-2.87118,-3.0982 ,-3.29462,-3.46995,-3.6309 ,-3.76107,-3.88989,-4.00288)
            if (k==2) B = (
                -3.10028,-3.42225,-3.64042,-3.84325,-4.00316,-4.13225,-4.25664,-4.36164,-4.49194,-4.59363 \
                -2.43198,-2.75066,-2.98703,-3.17889,-3.36513,-3.49435,-3.64109,-3.76414,-3.8762 ,-3.99358 \
                -2.07142,-2.39018,-2.63013,-2.83641,-3.01708,-3.16396,-3.30741,-3.43924,-3.55825,-3.68114)
            if (k==3) B = (
                -3.00871,-3.29276,-3.52116,-3.68947,-3.81867,-3.94979,-4.06436,-4.18565,-4.28861,-4.37955 \
                -2.36304,-2.63664,-2.86037,-3.03888,-3.16861,-3.31777,-3.4429 ,-3.54766,-3.65485,-3.76049 \
                -2.00831,-2.28933,-2.50595,-2.69041,-2.83741,-2.98112,-3.10822,-3.22499,-3.32096,-3.43545)
            if (k==4) B = (
                -3.00469,-3.26111,-3.45743,-3.60608,-3.76128,-3.86692,-3.94604,-4.05839,-4.16653,-4.24111 \
                -2.33238,-2.59496,-2.79945,-2.9735 ,-3.11097,-3.2225 ,-3.33309,-3.43588,-3.54444,-3.62254 \
                -1.97417,-2.24758,-2.45778,-2.62709,-2.77434,-2.89658,-3.0029 ,-3.12035,-3.21749,-3.30591)
            if (k==5) B = (
                -2.98685,-3.23316,-3.42661,-3.55741,-3.68509,-3.79058,-3.90019,-3.99035,-4.08989,-4.1701  \
                -2.31745,-2.57237,-2.77809,-2.92126,-3.05774,-3.18318,-3.27405,-3.37919,-3.47482,-3.56121 \
                -1.96448,-2.2257 ,-2.43366,-2.58093,-2.72721,-2.85165,-2.96099,-3.06389,-3.16189,-3.24555)
        }
        else { // n>=1000
            if (k==1) B = (
                -3.23764,-3.59173,-3.85336,-4.0583 ,-4.21101,-4.35501,-4.48991,-4.59941,-4.70568,-4.81266 \
                -2.58009,-2.9286 ,-3.20716,-3.42556,-3.61816,-3.76511,-3.93148,-4.04437,-4.16795,-4.27372 \
                -2.22325,-2.59103,-2.86825,-3.10536,-3.29575,-3.45961,-3.6288 ,-3.74991,-3.88434,-3.9968 )
            if (k==2) B = (
                -3.10513,-3.40641,-3.63596,-3.82943,-4.00141,-4.10874,-4.24518,-4.36941,-4.47675,-4.57808 \
                -2.41927,-2.74759,-2.97582,-3.17683,-3.35067,-3.5079 ,-3.64192,-3.76298,-3.88879,-3.9964  \
                -2.07046,-2.39083,-2.63585,-2.83003,-3.01836,-3.17335,-3.30697,-3.43725,-3.56603,-3.67783)
            if (k==3) B = (
                -3.03771,-3.3193 ,-3.50522,-3.6837 ,-3.81936,-3.94111,-4.07457,-4.16582,-4.27033,-4.36233 \
                -2.36914,-2.64359,-2.86459,-3.03056,-3.17932,-3.31424,-3.44624,-3.54739,-3.6552 ,-3.75265 \
                -2.00252,-2.29649,-2.51597,-2.68473,-2.84243,-2.97736,-3.11047,-3.21784,-3.33402,-3.42847)
            if (k==4) B = (
                -3.01191,-3.27634,-3.47498,-3.61148,-3.74011,-3.86386,-3.95855,-4.06205,-4.1629 ,-4.2342  \
                -2.3366 ,-2.60277,-2.80227,-2.97468,-3.10529,-3.23339,-3.33427,-3.43595,-3.53824,-3.61803 \
                -1.97311,-2.25349,-2.46026,-2.63308,-2.77237,-2.9027 ,-3.01563,-3.11598,-3.21908,-3.30184)
            if (k==5) B = (
                -2.99808,-3.26303,-3.42285,-3.56706,-3.69285,-3.80993,-3.91214,-3.97711,-4.06273,-4.15139 \
                -2.31276,-2.5889 ,-2.78662,-2.934  ,-3.05274,-3.18728,-3.2884 ,-3.37917,-3.46858,-3.55931 \
                -1.95523,-2.23706,-2.43915,-2.59319,-2.7281 ,-2.8627 ,-2.96611,-3.06685,-3.15656,-3.24319)
        }
    }
    return(B)
}

real rowvector __rals_fadf_cv(real scalar n, real scalar k,
                               real scalar model, real scalar r2)
{
    real matrix Bl, Bh
    real scalar nl, nh, w
    real rowvector lo, hi
    // bracket sample size
    if (n<=50) {
        nl = 50
        nh = 50
    }
    else if (n<=100) {
        nl = 50
        nh = 100
    }
    else if (n<=250) {
        nl = 100
        nh = 250
    }
    else if (n<=500) {
        nl = 250
        nh = 500
    }
    else if (n<=1000) {
        nl = 500
        nh = 1000
    }
    else {
        nl = 1000
        nh = 1000
    }
    Bl = __rals_fadf_block(nl, k, model)
    Bh = __rals_fadf_block(nh, k, model)
    // interpolate each row across r^2 grid
    lo = __rals_fadf_rinterp(Bl, r2)
    hi = __rals_fadf_rinterp(Bh, r2)
    if (nl==nh) return(lo)
    w = (n - nl) / (nh - nl)
    return((1-w)*lo + w*hi)
}

real rowvector __rals_fadf_rinterp(real matrix B, real scalar r2)
{
    // B is 3 x 10 ; grid is 0.1..1.0
    real scalar r210, ra, rb, wa
    real rowvector out
    if (r2<0.1) return(B[.,1]')
    r210 = r2*10
    if (r210>=10) return(B[.,10]')
    ra = floor(r210); rb = ceil(r210)
    if (ra==rb) return(B[.,ra]')
    wa = rb - r210
    out = wa*B[.,ra]' + (1-wa)*B[.,rb]'
    return(out)
}

//------------------------------------------------------------------------------
//  Long-run variance using Newey-West (Bartlett) -- needed for cointegration
//  rho^2 in the SNDE supplement.
//------------------------------------------------------------------------------
real matrix __rals_lrvar(real matrix u, real scalar bw)
{
    real scalar n, p, j, w
    real matrix omega
    n = rows(u); p = cols(u)
    if (bw<=0) bw = floor(4*(n/100)^(2/9))
    omega = quadcross(u,u)/n
    for (j=1; j<=bw; j++) {
        w = 1 - j/(bw+1)
        omega = omega + w*( quadcross(u[1::n-j,.],u[j+1::n,.])
                           +quadcross(u[j+1::n,.],u[1::n-j,.]) )/n
    }
    return(omega)
}

real scalar __rals_rho2_lr(real colvector res_aug, real colvector res_base)
{
    real matrix om
    real scalar r2
    om = __rals_lrvar((res_aug, res_base), 0)
    r2 = om[1,2]^2 / (om[1,1]*om[2,2])
    if (r2<0) r2 = 0
    if (r2>1) r2 = 1
    return(r2)
}

end

* End of top-level Mata block.  Library build itself is done by the
* program {bf:_rals_mata} defined at the top of this file.