*! ralsfkss 1.0.1  16may2026  Dr Merwan Roudane  <merwanroudane920@gmail.com>
*! RALS-Fourier KSS unit-root test
*  Reference:  Yilanci & Ozgur (2025) Politicka Ekonomie 73(3): 528-565.
*  Step 1: Fourier de-trending (Christopoulos & Leon-Ledesma 2010).
*  Step 2: KSS-type regression on residuals:  Δu_t = φ u^3_{t-1} + Σ ζ_i Δu_{t-i} + e_t
*  Step 3: Add RALS w-augmentation built from the stage-2 residuals.
*------------------------------------------------------------------------------

program define ralsfkss, rclass
    version 14.0
    capture mata: __rals_loaded()
    if _rc qui _rals_mata
    syntax varname(ts) [if] [in], [                ///
            TREND                                  ///
            MAXLags(integer 8)                     ///
            IC(string)                             ///
            FMAX(integer 5)                        ///
            FREQuency(real -1)                     ///
            Graph                                  ///
            noHEADer                               ]

    if `fmax' < 1 | `fmax' > 5 {
        di as error "fmax() must be between 1 and 5"
        exit 198
    }
    marksample touse
    qui count if `touse'
    if r(N) < 30 {
        di as error "sample too small for RALS-FKSS (N=" r(N) ")"
        exit 2001
    }
    capture tsset
    if _rc {
        di as error "data must be tsset before running ralsfkss"
        exit 459
    }
    local timevar `r(timevar)'

    if "`ic'"=="" local ic "aic"
    local ic = lower("`ic'")
    if !inlist("`ic'","aic","bic","tstat") {
        di as error "ic() must be aic, bic, or tstat"
        exit 198
    }
    local model = cond("`trend'"!="", 2, 1)

    tempname res cv
    mata: _rals_fkss_engine("`varlist'", "`touse'", `model', `maxlags', "`ic'", `fmax', `frequency', "`res'", "`cv'")

    local tauKSS  = `res'[1,1]
    local tauRALS = `res'[1,2]
    local rho2    = `res'[1,3]
    local T       = `res'[1,4]
    local lag     = `res'[1,5]
    local kfreq   = `res'[1,6]

    if "`header'"=="" {
        _rals_print_header "RALS-FKSS" "Yilanci & Ozgur (2025)" "`varlist'" `T' `lag' "`ic'" `model'
    }
    _rals_print_two_row "Fourier KSS (stage 1)" "RALS-FKSS (stage 2)"        ///
        `tauKSS' `tauRALS' `rho2'                                            ///
        `cv'[1,1] `cv'[1,2] `cv'[1,3] `cv'[2,1] `cv'[2,2] `cv'[2,3]
    di as text "  Optimal Fourier frequency k = " as result `kfreq'
    local concl_f = cond(`tauKSS'  < `cv'[1,2], "{bf:reject} H0 (FKSS, 5%)",            ///
                                                 "fail to reject H0 (FKSS, 5%)")
    local concl_r = cond(`tauRALS' < `cv'[2,2], "{bf:reject} H0 (RALS-FKSS, 5%)",       ///
                                                 "fail to reject H0 (RALS-FKSS, 5%)")
    di as text "  Decision (5%)"
    di as text "    {c -}{c -} FKSS      : " as result "`concl_f'"
    di as text "    {c -}{c -} RALS-FKSS : " as result "`concl_r'"
    di as text "{hline 80}"
    di as text "  Stage-1 CVs : Christopoulos & Leon-Ledesma (2010) Fourier-KSS table."
    di as text "  Stage-2 CVs : Yilanci & Ozgur (2025)."
    di as text "{hline 80}"

    return scalar tauKSS  = `tauKSS'
    return scalar tauRALS = `tauRALS'
    return scalar rho2    = `rho2'
    return scalar T       = `T'
    return scalar lag     = `lag'
    return scalar kfreq   = `kfreq'
    return scalar cv01_FKSS = `cv'[1,1]
    return scalar cv05_FKSS = `cv'[1,2]
    return scalar cv10_FKSS = `cv'[1,3]
    return scalar cv01      = `cv'[2,1]
    return scalar cv05      = `cv'[2,2]
    return scalar cv10      = `cv'[2,3]
    return local  test      = "RALS-FKSS"

    if "`graph'" != "" {
        _rals_graph_fadf "`varlist'" "`timevar'" "`touse'" `tauRALS' `cv'[2,2] `kfreq' `T'
    }
end

version 14.0
mata:
mata set matastrict off

void _rals_fkss_engine(string scalar yname, string scalar touse,
                       real scalar model, real scalar pmax, string scalar icmode,
                       real scalar fmax, real scalar fix_k,
                       string scalar rname, string scalar cvname)
{
    real colvector y, u, du, ulag, dep, s, c, t_, j
    real scalar T, k, kbest, ssrmin, popt, p, j_
    st_view(y, ., yname, touse)
    T = rows(y)
    real scalar kfrom, kto
    if (fix_k>0) {
        kfrom = fix_k
        kto   = fix_k
    }
    else {
        kfrom = 1
        kto   = fmax
    }

    // Stage 1: Fourier de-trending  (Christopoulos & Leon-Ledesma 2010)
    ssrmin = 1e+20; kbest = 1
    real colvector ubest
    ubest = J(T,1,.)
    for (k=kfrom; k<=kto; k++) {
        s = sin(2*pi()*(1::T)*k/T)
        c = cos(2*pi()*(1::T)*k/T)
        real matrix Z
        if (model==1) Z = J(T,1,1), s, c
        else          Z = J(T,1,1), (1::T), s, c
        real colvector bv, uv
        bv = invsym(quadcross(Z,Z))*quadcross(Z,y)
        uv = y - Z*bv
        real scalar ssrv
        ssrv = uv'uv
        if (ssrv < ssrmin) {
            ssrmin = ssrv
            kbest  = k
            ubest  = uv
        }
    }
    u = ubest

    // Stage 2: KSS regression  Δu_t = φ u^3_{t-1} + Σ ζ_i Δu_{t-i} + e_t
    du   = u[2::T] - u[1::T-1]
    ulag = u[1::T-1]
    real colvector u3
    u3 = ulag:^3
    // lag matrix of du
    real matrix lmat
    lmat = J(rows(du), pmax, 0)
    for (j_=1; j_<=pmax; j_++) lmat[j_+1::rows(du), j_] = du[1::rows(du)-j_]

    real matrix Xfull
    Xfull = u3, lmat
    real colvector depf
    depf = du
    if (pmax>0) {
        depf  = depf[pmax+1::rows(du)]
        Xfull = Xfull[pmax+1::rows(du), .]
    }

    // lag selection
    real colvector icvec, tvec
    real scalar bestic
    icvec = J(pmax+1,1,.); tvec = J(pmax+1,1,.)
    for (p=0; p<=pmax; p++) {
        real matrix Xp,XXi
        real colvector bp,ep,sep
        Xp  = Xfull[., 1::(1+p)]
        XXi = invsym(quadcross(Xp,Xp))
        bp  = XXi*quadcross(Xp,depf)
        ep  = depf - Xp*bp
        if (icmode=="bic") {
            icvec[p+1] = ln(ep'ep/rows(Xp)) + cols(Xp)*ln(rows(Xp))/rows(Xp)
        }
        else {
            icvec[p+1] = ln(ep'ep/rows(Xp)) + 2*cols(Xp)/rows(Xp)
        }
        if (p>0) {
            sep = sqrt(diagonal((ep'ep/(rows(Xp)-cols(Xp)))*XXi))
            tvec[p+1] = abs(bp[cols(Xp),1]/sep[cols(Xp),1])
        }
    }
    if (icmode=="tstat") {
        popt = 0
        for (p=pmax; p>=1; p--) {
            if (tvec[p+1] >= 1.645) {
                popt = p
                break
            }
        }
    }
    else {
        popt = 0
        bestic = icvec[1]
        for (p=1; p<=pmax; p++) {
            if (icvec[p+1] < bestic) {
                bestic = icvec[p+1]
                popt = p
            }
        }
    }

    real matrix Xs, XXi_s
    real colvector bs_, es_, ses_
    Xs    = Xfull[., 1::(1+popt)]
    XXi_s = invsym(quadcross(Xs,Xs))
    bs_   = XXi_s*quadcross(Xs,depf)
    es_   = depf - Xs*bs_
    real scalar sig2
    sig2  = (es_'es_)/(rows(Xs)-cols(Xs))
    ses_  = sqrt(diagonal(sig2*XXi_s))
    real scalar tauKSS
    tauKSS = bs_[1,1]/ses_[1,1]

    // RALS augmentation
    real scalar m2, m3, sig2A
    m2 = sum(es_:^2)/rows(es_)
    m3 = sum(es_:^3)/rows(es_)
    real matrix W, Xa, XaI
    real colvector ba, ea, sea
    W   = (es_:^2 :- m2), (es_:^3 :- m3 :- 3*m2*es_)
    Xa  = Xs, W
    XaI = invsym(quadcross(Xa,Xa))
    ba  = XaI*quadcross(Xa,depf)
    ea  = depf - Xa*ba
    sig2A = (ea'ea)/(rows(Xa)-cols(Xa))
    sea   = sqrt(diagonal(sig2A*XaI))
    real scalar tauRALS, rho2
    tauRALS = ba[1,1]/sea[1,1]
    rho2 = sig2A/sig2
    if (rho2<0) rho2=0
    if (rho2>1) rho2=1

    // CV.  Stage-2 RALS CVs use the RALS-LM table as a conservative proxy
    // (Yilanci-Ozgur 2025 §4.1 reports identical shape).
    real matrix crt, cvr
    crt = __rals_cv_lm()
    cvr = __rals_interp(crt, rho2)
    real rowvector cvKSS
    cvKSS = __rals_cv_fkss(T-popt, kbest, model)

    real matrix CV2
    CV2 = cvKSS \ cvr
    st_matrix(rname, (tauKSS, tauRALS, rho2, T, popt, kbest))
    st_matrix(cvname, CV2)
}
end