*! xtpanic 1.0.2  23jul2026   (xtflexur library)
*! PANIC panel unit root test (Bai & Ng 2004): tests the idiosyncratic
*! component after extracting common factors by principal components.
*!
*! Author : Dr Merwan Roudane  merwanroudane920@gmail.com
*!          https://github.com/merwanroudane
*! Part of the -xtflexur- library : factor-augmented panel time-series tests.
*!
*! Faithful port of the GAUSS routine BNG_PANIC (appl1_PANIC.gss / TSPDLIB) by
*!   Saban Nazlioglu.  Non-commercial use.
*!
*! References:
*!   Bai, J., Ng, S. (2004). A PANIC attack on unit roots and cointegration.
*!     Econometrica 72: 1127-1177. <doi:10.1111/j.1468-0262.2004.00528.x>
*!   Bai, J., Ng, S. (2002). Determining the number of factors in approximate
*!     factor models. Econometrica 70: 191-221.
*!
*! Step -> code map (see -help xtpanic methods-):
*!   number of common factors (Bai-Ng 2002 IC)  -> _xtp_fnumber
*!   principal-components factor extraction      -> _xtp_pca
*!   idiosyncratic ADF (no deterministics)       -> _xtp_adfnc
*!   pooled Fisher tests P and Pm                -> _xtp_run
*!   MC limiting distribution p-values           -> _xtp_pval

program xtpanic, rclass
    version 14.0

    syntax varname(numeric ts) [if] [in] [ ,   ///
        MODel(string)                           ///
        Kmax(integer 5)                         ///
        ICfactor(integer 2)                     ///
        Pmax(integer 3)                         ///
        IClag(integer 3)                        ///
    ]

    * deterministic model: 0 = constant, 1 = constant + trend
    if "`model'" == "" local model "constant"
    local model = lower("`model'")
    if      inlist("`model'","c","const","constant")  local p 0
    else if inlist("`model'","t","trend","ct")        local p 1
    else {
        di as error "xtpanic: model() must be constant or trend."
        exit 198
    }

    * panel structure
    quietly xtset
    local pv "`r(panelvar)'"
    local tv "`r(timevar)'"
    if "`pv'" == "" | "`tv'" == "" {
        di as error "xtpanic: data must be {help xtset}."
        exit 459
    }
    if "`r(balanced)'" != "strongly balanced" {
        di as error "xtpanic: requires a strongly balanced panel."
        exit 459
    }

    marksample touse
    markout `touse' `pv' `tv'
    quietly count if `touse'
    local NT = r(N)
    quietly levelsof `pv' if `touse', local(ids)
    local N : word count `ids'
    local T = `NT'/`N'
    if `T' != int(`T') | `N'*`T' != `NT' {
        di as error "xtpanic: unbalanced over the selected sample."
        exit 459
    }
    if `T' < 10 {
        di as error "xtpanic: need at least 10 time periods (found `T')."
        exit 2001
    }

    tempname R OUT
    mata: _xtp_run("`varlist'", "`touse'", `N', `T', `p', `kmax', `icfactor', `pmax', `iclag')
    matrix `R'   = r(pooled)
    matrix `OUT' = r(units)

    * ----- display -----------------------------------------------------------
    local mdesc = cond(`p'==0,"Constant","Constant and trend")
    local icfd = cond(`icfactor'==1,"PCp",cond(`icfactor'==2,"ICp","AIC/BIC"))
    di ""
    di as text "PANIC panel unit root test (Bai & Ng 2004)"
    di as text "{hline 60}"
    di as text "Panel var" _col(14) ": " as result "`pv'" ///
        as text _col(40) "N" _col(50) "= " as result %5.0f `N'
    di as text "Time var" _col(14) ": " as result "`tv'" ///
        as text _col(40) "T" _col(50) "= " as result %5.0f `T'
    di as text "Model" _col(14) ": " as result "`mdesc'" ///
        as text _col(40) "Factors" _col(50) "= " as result %3.0f `R'[3,1] ///
        as text " (`icfd')"
    di as text "H0: idiosyncratic component has a unit root (non-stationary)"
    di as text "{hline 60}"
    di as text "  Pooled test" _col(24) "Statistic" _col(40) "p-value"
    di as text "{hline 60}"
    local Pchi = `R'[1,1]
    local Pchip = `R'[1,2]
    local Pm   = `R'[2,1]
    local Pmp  = `R'[2,2]
    local s1 = cond(`Pchip'<0.01,"***",cond(`Pchip'<0.05,"**",cond(`Pchip'<0.10,"*","")))
    local s2 = cond(`Pmp'<0.01,"***",cond(`Pmp'<0.05,"**",cond(`Pmp'<0.10,"*","")))
    di as text "  P  (Fisher, chi2)" _col(24) as result %9.3f `Pchi' _col(40) %7.3f `Pchip' "  `s1'"
    di as text "  Pm (standardized)" _col(24) as result %9.3f `Pm' _col(40) %7.3f `Pmp' "  `s2'"
    di as text "{hline 60}"
    di as text "  Per-unit idiosyncratic ADF:"
    di as text "   id" _col(12) "ADF" _col(24) "p-value" _col(36) "lags"
    forvalues i = 1/`N' {
        di as text "  " %4.0f `OUT'[`i',1] _col(9) as result %9.3f `OUT'[`i',2] ///
            _col(24) %7.3f `OUT'[`i',3] _col(37) %3.0f `OUT'[`i',4]
    }
    di as text "{hline 60}"
    di as text "Note: reject H0 (=> stationary idiosyncratic component) for large"
    di as text "      P/Pm.  *** p<.01, ** p<.05, * p<.10."

    * ----- returns -----------------------------------------------------------
    return scalar P     = `Pchi'
    return scalar P_p   = `Pchip'
    return scalar Pm    = `Pm'
    return scalar Pm_p  = `Pmp'
    return scalar nf    = `R'[3,1]
    return scalar N     = `N'
    return scalar T     = `T'
    return matrix units = `OUT'
    return local  model "`mdesc'"
    return local  cmd   "xtpanic"
end


* =========================================================================
* Mata engine
* =========================================================================
version 14.0
mata:

// ---- first difference (pads a leading zero row) -----------------------
real matrix _xtp_diff(real matrix x)
{
    real scalar t
    t = rows(x)
    return((J(1,cols(x),0)) \ (x[2..t,.] - x[1..t-1,.]))
}

// ---- lag by 1 (pads a leading zero row) -------------------------------
real matrix _xtp_lagn(real matrix x)
{
    real scalar t
    t = rows(x)
    return((J(1,cols(x),0)) \ x[1..t-1,.])
}

// ---- principal-components factors ; returns Fhat, lambda, ehat --------
void _xtp_pca(real matrix X, real scalar nf, real matrix Fhat, real matrix lambda, real matrix ehat)
{
    real scalar t, n
    real matrix evec, A
    real rowvector eval
    real colvector ord
    t = rows(X)
    n = cols(X)
    if (t < n) {
        A = X*X'
        symeigensystem(A, evec, eval)
        ord = order(eval', -1)
        Fhat = sqrt(t) * evec[., ord[1..nf]]
        lambda = luinv(cross(Fhat,Fhat)) * cross(Fhat,X)
        ehat = X - Fhat*lambda
    }
    else {
        A = X'X
        symeigensystem(A, evec, eval)
        ord = order(eval', -1)
        lambda = sqrt(n) * evec[., ord[1..nf]]
        Fhat = X*lambda/n
        ehat = X - Fhat*lambda'
    }
}

// ---- number of factors (Bai-Ng 2002) ; returns (ICp1, ICp2) -----------
real rowvector _xtp_fnumber(real matrix X, real scalar kmax, real scalar ic)
{
    real scalar t, n, sigma2, k, mn, icp1, icp2
    real matrix Fhat, lambda, ehat, ICp
    real colvector V
    t = rows(X)
    n = cols(X)
    _xtp_pca(X, kmax, Fhat, lambda, ehat)
    sigma2 = sum(ehat:^2)/(n*t)
    mn = min((n,t))
    V = J(kmax,1,0)
    ICp = J(kmax,2,0)
    k = 1
    while (k <= kmax) {
        _xtp_pca(X, k, Fhat, lambda, ehat)
        V[k] = sum(ehat:^2)/(n*t)
        if (ic == 1) {
            ICp[k,1] = V[k] + k*sigma2*(n+t)/(n*t)*ln((n*t)/(n+t))
            ICp[k,2] = V[k] + k*sigma2*(n+t)/(n*t)*ln(mn)
        }
        else if (ic == 2) {
            ICp[k,1] = ln(V[k]) + k*(n+t)/(n*t)*ln((n*t)/(n+t))
            ICp[k,2] = ln(V[k]) + k*(n+t)/(n*t)*ln(mn)
        }
        else {
            ICp[k,1] = V[k] + k*sigma2*(2/t)
            ICp[k,2] = V[k] + k*sigma2*(ln(t)/t)
        }
        k = k + 1
    }
    icp1 = _xtp_argmin(ICp[,1])
    icp2 = _xtp_argmin(ICp[,2])
    return((icp1, icp2))
}

real scalar _xtp_argmin(real colvector v)
{
    real scalar i, mn, mi
    mn = v[1]
    mi = 1
    i = 2
    while (i <= rows(v)) {
        if (v[i] < mn) {
            mn = v[i]
            mi = i
        }
        i = i + 1
    }
    return(mi)
}

// ---- OLS: returns b, se ------------------------------------------------
void _xtp_ols(real colvector y, real matrix X, real colvector b, real colvector se)
{
    real matrix m
    real colvector e
    real scalar sig2
    m = luinv(cross(X,X))
    b = m*cross(X,y)
    e = y - X*b
    sig2 = cross(e,e)/(rows(y)-cols(X))
    se = sqrt(diagonal(m)*sig2)
}

// ---- ADF without deterministics ; returns (tau, lags) -----------------
real rowvector _xtp_adfnc(real colvector y, real scalar pmax, real scalar ic)
{
    real scalar t, p, j, tau, plag, n, k, LL
    real colvector dy, ly, dep, y1, b, se, taup, aicp, sicp, tstatp
    real matrix lmat, X
    t = rows(y)
    dy = _xtp_diff(y)
    ly = _xtp_lagn(y)
    taup = J(pmax+1,1,0)
    aicp = J(pmax+1,1,0)
    sicp = J(pmax+1,1,0)
    tstatp = J(pmax+1,1,0)
    p = 0
    while (p <= pmax) {
        lmat = J(rows(dy), p+1, 0)
        j = 1
        while (j <= p) {
            lmat[.,j] = _xtp_lagn_k(dy, j)
            j = j + 1
        }
        dep = dy[(p+2)..rows(dy)]
        y1  = ly[(p+2)..rows(ly)]
        if (p == 0) X = y1
        else        X = y1 , lmat[(p+2)..rows(lmat), 1..p]
        _xtp_ols(dep, X, b, se)
        taup[p+1] = b[1]/se[1]
        n = rows(X)
        k = cols(X)
        LL = -n/2*(1 + ln(2*pi()) + ln(cross(dep-X*b, dep-X*b)/n))
        aicp[p+1] = (2*k - 2*LL)/n
        sicp[p+1] = (k*ln(n) - 2*LL)/n
        tstatp[p+1] = abs(b[cols(X)]/se[cols(X)])
        p = p + 1
    }
    plag = _xtp_getlag(ic, pmax, aicp, sicp, tstatp)
    tau = taup[plag]
    return((tau, plag-1))
}

// lag a vector by k (pads leading zeros)
real colvector _xtp_lagn_k(real colvector x, real scalar k)
{
    real scalar t
    t = rows(x)
    return((J(k,1,0)) \ x[1..t-k])
}

// ---- lag selection -----------------------------------------------------
real scalar _xtp_getlag(real scalar ic, real scalar pmax, real colvector aicp, real colvector sicp, real colvector tstatp)
{
    real scalar p, j
    if (ic == 1) return(_xtp_argmin(aicp))
    if (ic == 2) return(_xtp_argmin(sicp))
    j = pmax + 1
    p = 1
    while (j >= 1) {
        if (abs(tstatp[j]) > 1.645) {
            p = j
            break
        }
        j = j - 1
    }
    if (j < 1) p = 1
    return(p)
}

// ---- ADF-no-constant p-value via the finite-sample response surface ---
//   CV(p,T) = c1 + c2/T + c3/T^2  (tspdlib panel.txt); nearest-CV lookup
real scalar _xtp_pval(real scalar stat, real scalar t)
{
    real matrix coe
    real colvector prob, par, qval
    real scalar j
    coe = _xtp_coeadfnc()
    prob = _xtp_dfprob()
    par = (1 \ 1/t \ 1/t^2)
    qval = coe*par
    if (stat < qval[1])   return(0.0001)
    if (stat > qval[399]) return(0.9999)
    j = _xtp_argmin(abs(qval :- stat))
    return(prob[j])
}

// ---- driver ------------------------------------------------------------
void _xtp_run(string scalar vname, string scalar touse, real scalar N, real scalar T, real scalar model, real scalar kmax, real scalar ick, real scalar pmax, real scalar icp)
{
    real colvector yv, adfe, lags, pval
    real matrix Y, x, Fhat, lambda, ehat, e, units, pooled
    real rowvector nf, ar
    real scalar Pchi, Pm, i, nfac

    yv = st_data(., vname, touse)
    // data are xtset-sorted (panel, time): row j*T is panel j
    Y = rowshape(yv, N)'          // T x N (column j = panel j)

    // difference and (for trend) demean
    x = _xtp_diff(Y)
    x = x[2..T, .]                // drop the padded first difference
    if (model == 1) x = x :- mean(x)

    nf = _xtp_fnumber(x, kmax, ick)
    nfac = nf[2]                   // ICp2
    _xtp_pca(x, nfac, Fhat, lambda, ehat)
    e = _xtp_cumsum(ehat)

    adfe = J(N,1,.)
    lags = J(N,1,.)
    pval = J(N,1,.)
    i = 1
    while (i <= N) {
        ar = _xtp_adfnc(e[.,i], pmax, icp)
        adfe[i] = ar[1]
        lags[i] = ar[2]
        // idiosyncratic ADF uses the DF-no-constant distribution for both
        // models (deterministics are removed by differencing/demeaning), as
        // in Bai-Ng (2004) / TSPDLIB and reproduced in the paper's Table 3.
        pval[i] = _xtp_pval(adfe[i], rows(e))
        i = i + 1
    }

    Pchi = -2*sum(ln(pval))
    Pm   = (Pchi - 2*N) / sqrt(4*N)

    units = (1::N) , adfe , pval , lags
    pooled = (Pchi , chi2tail(2*N, Pchi)) \ (Pm , 1-normal(Pm)) \ (nfac , .)
    st_matrix("r(units)", units)
    st_matrix("r(pooled)", pooled)
}

real matrix _xtp_cumsum(real matrix X)
{
    real matrix S
    real scalar j
    S = X
    j = 1
    while (j <= cols(X)) {
        S[.,j] = runningsum(X[.,j])
        j = j + 1
    }
    return(S)
}




// ---- ADF-no-constant response surface (tspdlib panel.txt) ----
// coe_adfnc[399,3] : CV(p,T) = c1 + c2/T + c3/T^2 ; df_prob[399] p-levels
real matrix _xtp_coeadfnc()
{
    string scalar s
    s = "-3.02200 -3.92700 -74.51000 -2.80300 -2.71300 -65.45000 -2.66700 -2.63600 -44.00000 -2.56800 -2.27500 -38.83000 -2.48800 -2.01900 -35.04000 -2.42200 -1.77100 -32.90000 -2.36500 -1.47700 -33.03000 -2.31500 -1.32400 -32.09000 -2.27000 -1.14100 -31.18000 -2.22900 -0.96320 -31.40000 -2.19200 -0.84700 -30.82000 -2.15700 -0.81540 -27.79000 -2.12400 -0.65930 -28.70000 -2.09400 -0.67410 -25.60000 -2.06500 -0.56160 -25.93000 -2.03700 -0.58010 -23.27000 -2.01200 -0.49750 -23.24000 -1.98700 -0.40110 -23.84000 -1.96400 -0.35030 -23.26000 -1.94200 -0.35700 -21.07000 -1.92000 -0.33320 -20.00000 -1.89900 -0.36500 -17.65000 -1.87900 -0.29940 -17.95000 -1.86000 -0.29020 -17.07000 -1.84100 -0.30300 -14.85000 -1.82400 -0.28420 -13.97000 -1.80600 -0.27880 -12.98000 -1.78900 -0.27220 -11.97000 -1.77300 -0.20610 -12.44000 -1.75700 -0.17650 -12.28000 -1.74100 -0.13890 -12.45000 -1.72600 -0.09571 -12.81000 -1.71100 -0.08605 -11.91000 -1.69700 -0.08214 -10.92000 -1.68300 -0.00952 -11.88000 -1.67000 0.01824 -11.70000 -1.65600 0.04615 -11.62000 -1.64300 0.09212 -11.88000 -1.63000 0.10860 -11.51000 -1.61800 0.16030 -11.93000 -1.60500 0.16000 -11.28000 -1.59300 0.17710 -11.14000 -1.58100 0.17950 -10.30000 -1.56900 0.17570 -9.43200 -1.55800 0.21280 -10.04000 -1.54600 0.21300 -9.43900 -1.53500 0.23560 -9.48300 -1.52400 0.23370 -8.46300 -1.51300 0.27900 -9.39100 -1.50200 0.26500 -8.52200 "
    s = s + "-1.49200 0.28960 -8.79000 -1.48100 0.27970 -8.12700 -1.47100 0.28590 -7.79100 -1.46000 0.28030 -6.95800 -1.45100 0.33010 -7.75100 -1.44100 0.35230 -7.73300 -1.43100 0.36230 -7.21700 -1.42200 0.36890 -6.67200 -1.41200 0.37660 -6.45600 -1.40300 0.37820 -6.16300 -1.39300 0.38470 -5.62100 -1.38400 0.39030 -5.31700 -1.37500 0.38020 -4.46200 -1.36600 0.40950 -4.78800 -1.35700 0.41270 -4.45900 -1.34900 0.43020 -4.38900 -1.34000 0.45890 -4.81500 -1.33200 0.49080 -5.11700 -1.32300 0.48630 -4.46300 -1.31500 0.51210 -4.70400 -1.30700 0.52690 -4.72700 -1.29800 0.54020 -4.74400 -1.29000 0.52540 -3.88300 -1.28200 0.50270 -2.98200 -1.27400 0.51920 -2.96300 -1.26600 0.52410 -2.52200 -1.25800 0.51010 -1.72600 -1.25000 0.52910 -1.86400 -1.24300 0.55360 -2.12600 -1.23500 0.57190 -2.29600 -1.22700 0.57390 -1.92600 -1.22000 0.58310 -1.65500 -1.21200 0.58250 -1.41300 -1.20500 0.60580 -1.49800 -1.19700 0.57680 -0.25140 -1.19000 0.56920 0.25050 -1.18200 0.55160 1.15500 -1.17500 0.54320 1.68300 -1.16800 0.54640 1.97900 -1.16100 0.57890 1.41100 -1.15400 0.58800 1.61300 -1.14700 0.61770 1.28100 -1.14000 0.60740 1.93000 -1.13300 0.62710 1.82400 -1.12600 0.63130 1.99400 -1.11900 0.65390 1.67000 -1.11200 0.67370 1.47700 -1.10600 0.67620 1.80200 -1.09900 0.67070 2.34300 -1.09200 0.68320 2.44100 "
    s = s + "-1.08600 0.68690 2.76800 -1.07900 0.69560 2.92900 -1.07200 0.68220 3.59800 -1.06600 0.68820 3.82400 -1.05900 0.69720 3.79800 -1.05300 0.71570 3.61800 -1.04600 0.71660 4.03600 -1.04000 0.72300 4.22100 -1.03300 0.72130 4.55700 -1.02700 0.70720 5.31300 -1.02000 0.69400 5.98300 -1.01400 0.70440 6.03100 -1.00800 0.71950 5.87100 -1.00200 0.73630 5.92300 -0.99540 0.75060 5.93800 -0.98920 0.76780 5.82800 -0.98300 0.77890 5.80800 -0.97700 0.80160 5.60900 -0.97090 0.81190 5.80900 -0.96460 0.79600 6.49900 -0.95850 0.81720 6.32400 -0.95250 0.83930 6.00100 -0.94630 0.83520 6.34200 -0.94040 0.85300 6.31100 -0.93430 0.86090 6.43100 -0.92840 0.86450 6.87900 -0.92240 0.87610 6.85700 -0.91650 0.90080 6.51900 -0.91050 0.90150 6.78500 -0.90460 0.91300 6.79000 -0.89870 0.92590 6.78900 -0.89280 0.92800 7.19000 -0.88690 0.93310 7.54600 -0.88100 0.93460 7.97400 -0.87510 0.93610 8.17800 -0.86930 0.95700 8.04200 -0.86360 0.97080 8.11900 -0.85780 0.98510 8.14900 -0.85200 0.99010 8.37700 -0.84610 0.99140 8.60400 -0.84040 1.01200 8.49900 -0.83470 1.03400 8.38600 -0.82890 1.04200 8.54400 -0.82300 1.03900 9.07900 -0.81740 1.06700 8.72400 -0.81180 1.10000 8.25100 -0.80610 1.11300 8.43300 -0.80040 1.12600 8.58200 -0.79460 1.12700 8.80400 -0.78910 1.17000 8.20200 "
    s = s + "-0.78350 1.20000 7.82100 -0.77790 1.21600 7.75600 -0.77220 1.22800 7.90300 -0.76660 1.23800 8.19900 -0.76100 1.26400 7.92500 -0.75540 1.29000 7.73900 -0.74980 1.31700 7.43100 -0.74420 1.34100 7.23400 -0.73850 1.37100 6.84900 -0.73290 1.39500 6.62300 -0.72730 1.43200 5.93800 -0.72180 1.46200 5.71000 -0.71620 1.47600 5.84300 -0.71060 1.50000 5.58400 -0.70500 1.52900 5.21000 -0.69940 1.54300 5.42200 -0.69360 1.55700 5.42400 -0.68810 1.59400 4.96900 -0.68250 1.61700 4.89600 -0.67690 1.64300 4.73700 -0.67140 1.66900 4.74000 -0.66580 1.71000 4.12700 -0.66010 1.73000 4.20900 -0.65440 1.75700 3.98600 -0.64880 1.78600 3.85900 -0.64310 1.80100 3.93700 -0.63750 1.85800 2.81900 -0.63180 1.88500 2.59600 -0.62600 1.89300 2.88000 -0.62010 1.89600 3.31700 -0.61440 1.93200 2.79100 -0.60870 1.96000 2.60000 -0.60300 1.99400 2.12200 -0.59720 2.00600 2.20300 -0.59140 2.03700 1.78600 -0.58570 2.08600 0.95880 -0.57980 2.10000 0.91940 -0.57380 2.11000 1.08900 -0.56780 2.12500 1.06100 -0.56210 2.15900 0.82260 -0.55620 2.19400 0.21650 -0.55030 2.21100 0.21700 -0.54440 2.23800 0.12570 -0.53870 2.28500 -0.51150 -0.53280 2.31900 -0.83040 -0.52690 2.34800 -1.12500 -0.52100 2.38900 -1.90200 -0.51490 2.41300 -2.14800 -0.50890 2.44000 -2.39400 -0.50280 2.45500 -2.29100 "
    s = s + "-0.49680 2.49300 -2.79900 -0.49060 2.49300 -2.46600 -0.48460 2.53100 -2.98000 -0.47850 2.56700 -3.40500 -0.47240 2.60300 -4.03500 -0.46620 2.63800 -4.65600 -0.46010 2.67000 -5.08400 -0.45400 2.70400 -5.66300 -0.44770 2.72700 -5.98600 -0.44150 2.74700 -6.18900 -0.43540 2.78900 -6.86100 -0.42910 2.82300 -7.64400 -0.42280 2.84900 -8.07400 -0.41650 2.88200 -8.62900 -0.41020 2.89500 -8.74900 -0.40370 2.89400 -8.35800 -0.39720 2.90200 -8.28600 -0.39080 2.92200 -8.56600 -0.38430 2.93600 -8.84700 -0.37770 2.95200 -9.09400 -0.37120 2.96700 -9.26600 -0.36480 2.99800 -9.93600 -0.35810 3.00600 -10.01000 -0.35160 3.02000 -10.15000 -0.34490 3.04300 -10.87000 -0.33830 3.06500 -11.30000 -0.33190 3.11200 -12.52000 -0.32520 3.11300 -12.34000 -0.31850 3.13300 -12.76000 -0.31190 3.15700 -13.21000 -0.30530 3.17500 -13.58000 -0.29860 3.18800 -13.72000 -0.29190 3.21100 -14.05000 -0.28510 3.22800 -14.35000 -0.27840 3.25100 -14.85000 -0.27140 3.23200 -14.27000 -0.26460 3.24400 -14.41000 -0.25780 3.26300 -14.94000 -0.25090 3.27700 -15.34000 -0.24380 3.25200 -14.62000 -0.23710 3.27600 -15.07000 -0.23020 3.29200 -15.31000 -0.22310 3.27600 -14.92000 -0.21620 3.29600 -15.56000 -0.20920 3.28800 -15.17000 -0.20230 3.31100 -15.62000 -0.19540 3.33900 -16.25000 -0.18820 3.31700 -15.51000 -0.18110 3.30200 -15.03000 -0.17400 3.30800 -15.09000 "
    s = s + "-0.16690 3.31100 -14.95000 -0.15970 3.31200 -14.96000 -0.15250 3.32300 -15.23000 -0.14550 3.33000 -15.19000 -0.13830 3.34800 -15.64000 -0.13120 3.36800 -16.15000 -0.12400 3.36800 -16.05000 -0.11680 3.38900 -16.60000 -0.10930 3.35800 -15.80000 -0.10190 3.35000 -15.52000 -0.09462 3.35200 -15.53000 -0.08714 3.35000 -15.47000 -0.07975 3.34500 -15.32000 -0.07223 3.35000 -15.59000 -0.06471 3.34500 -15.43000 -0.05731 3.35600 -15.71000 -0.04975 3.34600 -15.37000 -0.04210 3.34600 -15.45000 -0.03445 3.33300 -15.17000 -0.02671 3.31600 -14.82000 -0.01906 3.32500 -15.20000 -0.01125 3.30100 -14.59000 -0.00369 3.31800 -15.05000 0.00386 3.32400 -15.21000 0.01152 3.32200 -15.11000 0.01924 3.31500 -14.81000 0.02709 3.31300 -14.81000 0.03492 3.31000 -14.82000 0.04276 3.31400 -15.00000 0.05045 3.34300 -15.70000 0.05831 3.35700 -16.06000 0.06624 3.36400 -16.22000 0.07408 3.38800 -16.83000 0.08201 3.39000 -16.84000 0.08999 3.38600 -16.62000 0.09795 3.40400 -17.01000 0.10610 3.40500 -17.09000 0.11410 3.42100 -17.47000 0.12210 3.43900 -18.00000 0.13040 3.43800 -18.02000 0.13850 3.47000 -18.80000 0.14680 3.47300 -18.82000 0.15540 3.44100 -17.93000 0.16390 3.42000 -17.43000 0.17220 3.43200 -17.85000 0.18050 3.44800 -18.20000 0.18910 3.43100 -17.72000 0.19770 3.43400 -17.75000 0.20620 3.42800 -17.68000 0.21480 3.43400 -17.78000 "
    s = s + "0.22350 3.45100 -18.32000 0.23240 3.42800 -17.79000 0.24120 3.42100 -17.66000 0.25020 3.41000 -17.57000 0.25900 3.42000 -17.94000 0.26790 3.42600 -18.08000 0.27710 3.40000 -17.48000 0.28600 3.42700 -18.23000 0.29510 3.43700 -18.64000 0.30430 3.41700 -17.97000 0.31340 3.43700 -18.30000 0.32250 3.45500 -18.86000 0.33210 3.41600 -18.02000 0.34170 3.40300 -17.78000 0.35120 3.39700 -17.54000 0.36090 3.38800 -17.49000 0.37050 3.38300 -17.28000 0.38010 3.38400 -17.27000 0.38990 3.38100 -17.33000 0.39960 3.41500 -18.25000 0.40940 3.42700 -18.79000 0.41920 3.44900 -19.15000 0.42930 3.43900 -18.88000 0.43940 3.44000 -18.98000 0.44940 3.45500 -19.42000 0.45980 3.46200 -19.65000 0.47010 3.46500 -19.55000 0.48060 3.45600 -19.26000 0.49090 3.46600 -19.22000 0.50160 3.47000 -19.55000 0.51240 3.48200 -19.90000 0.52300 3.50500 -20.30000 0.53400 3.49500 -20.04000 0.54520 3.48400 -20.11000 0.55630 3.48400 -19.78000 0.56740 3.48900 -19.52000 0.57870 3.48300 -19.11000 0.59020 3.51100 -20.12000 0.60160 3.54200 -20.98000 0.61370 3.50700 -20.14000 0.62540 3.52700 -20.67000 0.63750 3.51700 -20.21000 0.64950 3.52500 -20.40000 0.66190 3.53200 -20.69000 0.67410 3.55200 -20.89000 0.68680 3.56700 -21.27000 0.69940 3.58300 -21.58000 0.71220 3.60500 -22.09000 0.72510 3.63900 -23.02000 0.73840 3.65100 -23.31000 "
    s = s + "0.75230 3.59400 -21.78000 0.76610 3.55800 -20.66000 0.78000 3.55700 -20.34000 0.79410 3.55400 -20.19000 0.80830 3.56400 -20.45000 0.82310 3.51000 -18.95000 0.83750 3.54700 -19.34000 0.85240 3.55500 -19.33000 0.86740 3.56000 -19.29000 0.88310 3.52100 -18.16000 0.89880 3.54600 -18.77000 0.91490 3.51900 -17.56000 0.93140 3.51900 -17.68000 0.94810 3.50700 -16.79000 0.96490 3.55500 -17.42000 0.98220 3.58400 -17.84000 1.00000 3.58700 -18.10000 1.01800 3.60700 -18.34000 1.03700 3.63400 -18.56000 1.05600 3.64000 -18.82000 1.07500 3.67000 -18.81000 1.09500 3.64400 -17.48000 1.11600 3.63000 -16.65000 1.13700 3.66900 -17.28000 1.15900 3.64900 -16.32000 1.18200 3.68800 -16.85000 1.20500 3.67200 -15.87000 1.22900 3.76100 -17.59000 1.25400 3.77000 -17.16000 1.28000 3.77500 -16.67000 1.30600 3.86000 -17.65000 1.33400 3.97300 -19.73000 1.36200 4.05300 -20.09000 1.39200 4.13200 -20.97000 1.42500 4.16300 -20.39000 1.45900 4.24600 -21.49000 1.49500 4.30200 -22.44000 1.53300 4.37300 -22.32000 1.57400 4.44800 -22.38000 1.61900 4.57000 -23.34000 1.66700 4.60900 -22.29000 1.72000 4.61800 -19.71000 1.77900 4.70200 -18.44000 1.84500 4.76700 -15.42000 1.92000 4.95500 -12.53000 2.01000 5.31900 -13.95000 2.12100 5.80200 -14.49000 2.27100 6.76700 -22.19000 2.51400 7.77700 -16.58000 "
    return(rowshape(strtoreal(tokens(s)), 399))
}

real colvector _xtp_dfprob()
{
    string scalar s
    s = "0.00250 0.00500 0.00750 0.01000 0.01250 0.01500 0.01750 0.02000 0.02250 0.02500 0.02750 0.03000 0.03250 0.03500 0.03750 0.04000 0.04250 0.04500 0.04750 0.05000 0.05250 0.05500 0.05750 0.06000 0.06250 0.06500 0.06750 0.07000 0.07250 0.07500 0.07750 0.08000 0.08250 0.08500 0.08750 0.09000 0.09250 0.09500 0.09750 0.10000 0.10250 0.10500 0.10750 0.11000 0.11250 0.11500 0.11750 0.12000 0.12250 0.12500 0.12750 0.13000 0.13250 0.13500 0.13750 0.14000 0.14250 0.14500 0.14750 0.15000 0.15250 0.15500 0.15750 0.16000 0.16250 0.16500 0.16750 0.17000 0.17250 0.17500 0.17750 0.18000 0.18250 0.18500 0.18750 0.19000 0.19250 0.19500 0.19750 0.20000 0.20250 0.20500 0.20750 0.21000 0.21250 0.21500 0.21750 0.22000 0.22250 0.22500 0.22750 0.23000 0.23250 0.23500 0.23750 0.24000 0.24250 0.24500 0.24750 0.25000 0.25250 0.25500 0.25750 0.26000 0.26250 0.26500 0.26750 0.27000 0.27250 0.27500 0.27750 0.28000 0.28250 0.28500 0.28750 0.29000 0.29250 0.29500 0.29750 0.30000 0.30250 0.30500 0.30750 0.31000 0.31250 0.31500 0.31750 0.32000 0.32250 0.32500 0.32750 0.33000 0.33250 0.33500 0.33750 0.34000 0.34250 0.34500 0.34750 0.35000 0.35250 0.35500 0.35750 0.36000 0.36250 0.36500 0.36750 0.37000 0.37250 0.37500 "
    s = s + "0.37750 0.38000 0.38250 0.38500 0.38750 0.39000 0.39250 0.39500 0.39750 0.40000 0.40250 0.40500 0.40750 0.41000 0.41250 0.41500 0.41750 0.42000 0.42250 0.42500 0.42750 0.43000 0.43250 0.43500 0.43750 0.44000 0.44250 0.44500 0.44750 0.45000 0.45250 0.45500 0.45750 0.46000 0.46250 0.46500 0.46750 0.47000 0.47250 0.47500 0.47750 0.48000 0.48250 0.48500 0.48750 0.49000 0.49250 0.49500 0.49750 0.50000 0.50250 0.50500 0.50750 0.51000 0.51250 0.51500 0.51750 0.52000 0.52250 0.52500 0.52750 0.53000 0.53250 0.53500 0.53750 0.54000 0.54250 0.54500 0.54750 0.55000 0.55250 0.55500 0.55750 0.56000 0.56250 0.56500 0.56750 0.57000 0.57250 0.57500 0.57750 0.58000 0.58250 0.58500 0.58750 0.59000 0.59250 0.59500 0.59750 0.60000 0.60250 0.60500 0.60750 0.61000 0.61250 0.61500 0.61750 0.62000 0.62250 0.62500 0.62750 0.63000 0.63250 0.63500 0.63750 0.64000 0.64250 0.64500 0.64750 0.65000 0.65250 0.65500 0.65750 0.66000 0.66250 0.66500 0.66750 0.67000 0.67250 0.67500 0.67750 0.68000 0.68250 0.68500 0.68750 0.69000 0.69250 0.69500 0.69750 0.70000 0.70250 0.70500 0.70750 0.71000 0.71250 0.71500 0.71750 0.72000 0.72250 0.72500 0.72750 0.73000 0.73250 0.73500 0.73750 0.74000 0.74250 0.74500 0.74750 0.75000 "
    s = s + "0.75250 0.75500 0.75750 0.76000 0.76250 0.76500 0.76750 0.77000 0.77250 0.77500 0.77750 0.78000 0.78250 0.78500 0.78750 0.79000 0.79250 0.79500 0.79750 0.80000 0.80250 0.80500 0.80750 0.81000 0.81250 0.81500 0.81750 0.82000 0.82250 0.82500 0.82750 0.83000 0.83250 0.83500 0.83750 0.84000 0.84250 0.84500 0.84750 0.85000 0.85250 0.85500 0.85750 0.86000 0.86250 0.86500 0.86750 0.87000 0.87250 0.87500 0.87750 0.88000 0.88250 0.88500 0.88750 0.89000 0.89250 0.89500 0.89750 0.90000 0.90250 0.90500 0.90750 0.91000 0.91250 0.91500 0.91750 0.92000 0.92250 0.92500 0.92750 0.93000 0.93250 0.93500 0.93750 0.94000 0.94250 0.94500 0.94750 0.95000 0.95250 0.95500 0.95750 0.96000 0.96250 0.96500 0.96750 0.97000 0.97250 0.97500 0.97750 0.98000 0.98250 0.98500 0.98750 0.99000 0.99250 0.99500 0.99750 "
    return(strtoreal(tokens(s))')
}

end