*! xtdhcoint v1.0.0 - Durbin-Hausman Panel Cointegration Tests *! Faithful Stata implementation of the original GAUSS code by Westerlund *! Date: February 2026 *! Author: Dr. Merwan Roudane (merwanroudane920@gmail.com) *! Reference: Westerlund, J. (2008). Panel cointegration tests of the Fisher effect. *! Journal of Applied Econometrics, 23: 193-233. /* ** VERSION 2.0.0 - COMPLETE REWRITE ** ** Changes from v1.0.x: ** - Faithful line-by-line replication of the original GAUSS code (dhh.txt) ** - Zero-padding before differencing (matches GAUSS cum() lines 19-22) ** - Robust balanced-panel extraction ** - Replaced 'mata clear' with targeted function drops ** - Removed bootstrap (not in original GAUSS code) ** ** GAUSS CODE REFERENCE (dhh.txt by Joakim Westerlund): ** cum() -> lines 10-38 : First-diff, defactor, cumulate ** gdh() -> lines 40-63 : Individual DH statistic ** pdh() -> lines 65-102 : Panel DH statistic ** fact() -> lines 104-158 : Factor number selection ** prin() -> lines 160-178 : Principal components ** gdh_panel() -> lines 180-207 : Main driver, standardization ** ** PAPER REFERENCE (Westerlund 2008, p.203-204): ** DHg = sum_i Si * (phi_tilde_i - phi_hat_i)^2 * sum_t(e_{it-1}^2) ** DHp = Sn * (phi_tilde - phi_hat)^2 * sum_i sum_t(e_{it-1}^2) ** Standardization via Theorem 1 moments: ** E(Bi)=5.5464, var(Bi)=36.7673, E(Ci)=0.5005, var(Ci)=0.3348 */ program define xtdhcoint, rclass sortpreserve version 14.0 syntax varlist(min=2 numeric ts) [if] [in], /// [Kmax(integer 5)] /// Maximum number of factors [Criterion(string)] /// Information criterion (ic/pc/aic/bic) [Penalty(integer 1)] /// Penalty type (1=default, 2, 3) [Bandwidth(integer -1)] /// Kernel bandwidth (-1=automatic) [PREdet(real 0)] /// Predetermined coefficient (0=estimate) [NOConstant] /// Suppress constant term [Level(cilevel)] /// Confidence level [NOTable] /// Suppress output table *-------------------------------------------------------------------------- * 1. Parse Variables and Check Panel Structure *-------------------------------------------------------------------------- marksample touse * Get variable names gettoken depvar indepvars : varlist local nvars : word count `indepvars' * Check panel structure qui xtset local panelvar `r(panelvar)' local timevar `r(timevar)' if "`panelvar'" == "" { di as error "Panel variable not set. Use xtset panelvar timevar first." exit 459 } *-------------------------------------------------------------------------- * 2. Set Default Options *-------------------------------------------------------------------------- * Information criterion if "`criterion'" == "" local criterion "ic" if !inlist("`criterion'", "ic", "pc", "aic", "bic") { di as error "criterion() must be one of: ic, pc, aic, bic" exit 198 } * Map criterion to code (matches GAUSS fact() parameter c) if "`criterion'" == "pc" local cri = 1 if "`criterion'" == "ic" local cri = 2 if "`criterion'" == "aic" local cri = 3 if "`criterion'" == "bic" local cri = 3 * BIC uses penalty type 2 if "`criterion'" == "bic" local penalty = 2 *-------------------------------------------------------------------------- * 3. Extract Balanced Panel Data *-------------------------------------------------------------------------- preserve qui keep if `touse' * Sort by panel and time sort `panelvar' `timevar' * Get panel IDs tempvar panel_id qui egen `panel_id' = group(`panelvar') qui sum `panel_id', meanonly local N = r(max) if `N' < 2 { di as error "At least 2 cross-sectional units required." restore exit 198 } * Get time periods per panel and enforce balance tempvar time_count qui bysort `panel_id': gen `time_count' = _N qui sum `time_count', meanonly local Tmin = r(min) local Tmax = r(max) if `Tmin' != `Tmax' { di as text "Note: Unbalanced panel. Truncating each unit to T = " `Tmin' } * Keep only first Tmin observations per unit for balance tempvar obs_within qui bysort `panel_id' (`timevar'): gen `obs_within' = _n qui keep if `obs_within' <= `Tmin' local T = `Tmin' if `T' < 10 { di as error "At least 10 time periods required." restore exit 198 } * Automatic bandwidth (Newey-West) - matches paper Section 3.4 if `bandwidth' < 0 { local bandwidth = floor(4 * (`T'/100)^(2/9)) } if `bandwidth' < 1 { local bandwidth = 1 } * Maximum factors validation if `kmax' < 1 { di as error "kmax() must be at least 1" restore exit 198 } if `kmax' > min(`N', `T') { local kmax = min(`N', `T') - 1 di as text "Note: kmax adjusted to " `kmax' } * Re-sort to ensure correct order for reshaping sort `panel_id' `timevar' *-------------------------------------------------------------------------- * 4. Verify data dimensions are valid for Mata *-------------------------------------------------------------------------- qui count local n_obs = r(N) local expected_obs = `N' * `T' if `n_obs' != `expected_obs' { di as error "Data dimension mismatch: expected " `expected_obs' " obs, found " `n_obs' restore exit 198 } *-------------------------------------------------------------------------- * 5. Pass to Mata for Computation *-------------------------------------------------------------------------- mata: _xtdh_main("`depvar'", "`indepvars'", "`panel_id'", /// `N', `T', `kmax', `penalty', `cri', /// `bandwidth', `predet') *-------------------------------------------------------------------------- * 6. Retrieve Results *-------------------------------------------------------------------------- tempname dhg dhp dhg_z dhp_z dhg_p dhp_p nf scalar `dhg' = r(dhg) scalar `dhp' = r(dhp) scalar `dhg_z' = r(dhg_z) scalar `dhp_z' = r(dhp_z) scalar `dhg_p' = r(dhg_p) scalar `dhp_p' = r(dhp_p) scalar `nf' = r(nf) restore *-------------------------------------------------------------------------- * 7. Display Results *-------------------------------------------------------------------------- if "`notable'" == "" { di "" di as text "{hline 78}" di as text "{bf:Westerlund (2008) Durbin-Hausman Panel Cointegration Tests}" di as text "{hline 78}" di "" di as text "H0: No cointegration" di as text "H1(DHg): Cointegration for at least some units" di as text "H1(DHp): Cointegration for all units (common AR parameter)" di "" di as text "{hline 78}" di as text "Panel dimensions:" _col(25) "N = " as result %5.0f `N' /// as text _col(45) "T = " as result %5.0f `T' di as text "Number of regressors:" _col(25) "K = " as result %5.0f `nvars' di as text "Estimated factors:" _col(25) "F = " as result %5.0f scalar(`nf') di as text "Kernel bandwidth:" _col(25) "M = " as result %5.0f `bandwidth' if `predet' != 0 { di as text "Predetermined coeff.:" _col(25) as result %7.4f `predet' } di as text "{hline 78}" di "" if `nvars' > 1 { di as error "Warning: Simulated moments are derived for K=1 (one regressor)." di as error " With K=`nvars', z-values may not follow N(0,1) under H0." di as error " Consider running separate bivariate tests." di "" } * Compute averages for display local dhg_avg = scalar(`dhg') / `N' local dhp_avg = scalar(`dhp') / `N' * Theoretical means under H0 (Westerlund 2008, Remark 2, p.204) local mu_g = 5.5464 local mu_p_inv = 1 / 0.5005 // E(1/Ci) = 1/E(Ci) used in DHp * Results table di as text "{hline 78}" di as text _col(5) "Test" _col(22) "Avg.Stat" _col(35) "E(.) H0" /// _col(48) "Z-value" _col(61) "P-value" _col(73) "" di as text "{hline 78}" * DHg result local dhg_dec = cond(scalar(`dhg_p') < (100 - `level')/100, "Reject", "Fail to reject") di as text _col(5) "DHg" /// _col(20) as result %10.4f `dhg_avg' /// _col(33) as result %10.4f `mu_g' /// _col(46) as result %10.4f scalar(`dhg_z') /// _col(59) as result %10.4f scalar(`dhg_p') /// _col(73) as text "`dhg_dec'" * DHp result local dhp_dec = cond(scalar(`dhp_p') < (100 - `level')/100, "Reject", "Fail to reject") di as text _col(5) "DHp" /// _col(20) as result %10.4f `dhp_avg' /// _col(33) as result %10.4f `mu_p_inv' /// _col(46) as result %10.4f scalar(`dhp_z') /// _col(59) as result %10.4f scalar(`dhp_p') /// _col(73) as text "`dhp_dec'" di as text "{hline 78}" di "" di as text "Avg.Stat = average per-unit statistic (raw sum / N)." di as text "E(.) H0 = expected value under H0. Z = sqrt(N)*(Avg.Stat - E(.))/se." di "" * ---------------------------------------------------------------- * Decision Table: Z-value vs Critical Values * ---------------------------------------------------------------- * Critical values (right tail of standard normal) local cv10 = invnormal(0.90) // 1.2816 local cv05 = invnormal(0.95) // 1.6449 local cv01 = invnormal(0.99) // 2.3263 * Decision functions local dhg_d10 = cond(scalar(`dhg_z') > `cv10', "Reject", " -- ") local dhg_d05 = cond(scalar(`dhg_z') > `cv05', "Reject", " -- ") local dhg_d01 = cond(scalar(`dhg_z') > `cv01', "Reject", " -- ") local dhp_d10 = cond(scalar(`dhp_z') > `cv10', "Reject", " -- ") local dhp_d05 = cond(scalar(`dhp_z') > `cv05', "Reject", " -- ") local dhp_d01 = cond(scalar(`dhp_z') > `cv01', "Reject", " -- ") di as text "{bf:Decision: Z-value vs Critical Values (right tail)}" di as text "{hline 78}" di as text _col(5) "Test" _col(20) "Z-value" /// _col(35) "CV 10%" _col(48) "CV 5%" _col(60) "CV 1%" _col(73) "" di as text _col(35) as result %7.4f `cv10' /// _col(47) as result %7.4f `cv05' /// _col(59) as result %7.4f `cv01' di as text "{hline 78}" di as text _col(5) "DHg" /// _col(18) as result %10.4f scalar(`dhg_z') /// _col(35) as text "`dhg_d10'" /// _col(48) as text "`dhg_d05'" /// _col(60) as text "`dhg_d01'" di as text _col(5) "DHp" /// _col(18) as result %10.4f scalar(`dhp_z') /// _col(35) as text "`dhp_d10'" /// _col(48) as text "`dhp_d05'" /// _col(60) as text "`dhp_d01'" di as text "{hline 78}" di "" di as text "Reject = Z-value > Critical Value (evidence of cointegration)." di as text " -- = Fail to reject H0 (no evidence of cointegration)." di "" di as text "Reference: Westerlund (2008), J. of Applied Econometrics, 23: 193-233" di "" } *-------------------------------------------------------------------------- * 8. Store Results *-------------------------------------------------------------------------- return scalar N = `N' return scalar T = `T' return scalar K = `nvars' return scalar nf = scalar(`nf') return scalar bandwidth = `bandwidth' return scalar kmax = `kmax' return scalar dhg = scalar(`dhg') return scalar dhp = scalar(`dhp') return scalar dhg_z = scalar(`dhg_z') return scalar dhp_z = scalar(`dhp_z') return scalar dhg_p = scalar(`dhg_p') return scalar dhp_p = scalar(`dhp_p') return local depvar "`depvar'" return local indepvars "`indepvars'" return local criterion "`criterion'" end *============================================================================== * MATA FUNCTIONS * Line-by-line replication of GAUSS code (dhh.txt by Joakim Westerlund) *============================================================================== * Safely drop existing functions before redefining capture mata: mata drop _xtdh_main() capture mata: mata drop _xtdh_cum() capture mata: mata drop _xtdh_gdh() capture mata: mata drop _xtdh_pdh() capture mata: mata drop _xtdh_fact() capture mata: mata drop _xtdh_prin() capture mata: mata drop _xtdh_fejer() mata: mata set matastrict on /* ** Main driver function ** Replicates GAUSS gdh_panel() lines 180-207 */ void _xtdh_main(string scalar depvar, string scalar indepvars, string scalar panelid, real scalar N, real scalar T, real scalar kmax, real scalar pen, real scalar cri, real scalar bandwidth, real scalar predet) { real matrix e real scalar i, nf, dhg, dhg_raw, dhp_raw, dhg_z, dhp_z, dhg_p, dhp_p real scalar n_test // Simulated moments from Westerlund (2008) Remark 2, p.204 // Bi = (integral W(r)^2 dr)^{-1} // Ci = 1/Bi = integral W(r)^2 dr real scalar mu_g, var_g, mu_p, var_p mu_g = 5.5464 // E(Bi) var_g = 36.7673 // var(Bi) mu_p = 0.5005 // E(Ci) var_p = 0.3348 // var(Ci) // Step 1: Compute cumulated defactored residuals // Replicates GAUSS: {e, nf} = cum(x, y, kmax, pen, cri) _xtdh_cum(depvar, indepvars, panelid, N, T, kmax, pen, cri, predet, e, nf) // n_test = number of units to test (= N in general use) // GAUSS line 186: n = cols(y) n_test = cols(e) // Step 2: Compute DHg = sum of individual DH statistics // Replicates GAUSS gdh_panel() lines 190-195 dhg_raw = 0 for (i = 1; i <= n_test; i++) { dhg_raw = dhg_raw + _xtdh_gdh(e[., i], bandwidth) } // Step 3: Compute DHp = panel DH statistic dhp_raw = _xtdh_pdh(e, bandwidth) // Step 4: Standardize - GAUSS gdh_panel() lines 198-204 // GAUSS line 200: dhg = sqrt(n)*(dhg/n - mu_g)/sqrt(var_g) dhg_z = sqrt(n_test) * (dhg_raw / n_test - mu_g) / sqrt(var_g) // GAUSS line 204: dhp = sqrt(n)*(pdh(e,p)/n - mu_p^(-1))/sqrt(mu_p^(-4)*var_p) dhp_z = sqrt(n_test) * (dhp_raw / n_test - mu_p^(-1)) / sqrt(mu_p^(-4) * var_p) // P-values (right tail of standard normal) // Remark 2: "the computed value should be compared to the right tail" dhg_p = 1 - normal(dhg_z) dhp_p = 1 - normal(dhp_z) // Store results st_numscalar("r(dhg)", dhg_raw) st_numscalar("r(dhp)", dhp_raw) st_numscalar("r(dhg_z)", dhg_z) st_numscalar("r(dhp_z)", dhp_z) st_numscalar("r(dhg_p)", dhg_p) st_numscalar("r(dhp_p)", dhp_p) st_numscalar("r(nf)", nf) } /* ** Cumulated defactored residuals ** Replicates GAUSS cum() lines 10-38 ** ** Key: GAUSS prepends zeros to y and x before differencing (lines 19-20) ** This gives T first-differenced observations (not T-1) */ void _xtdh_cum(string scalar depvar, string scalar indepvars, string scalar panelid, real scalar N, real scalar T, real scalar kmax, real scalar pen, real scalar cri, real scalar predet, real matrix e, real scalar nf) { real matrix Y_raw, X_raw, Y_wide, de, F, Lambda real colvector panel real scalar i, K // Read data from Stata (sorted by panel then time) Y_raw = st_data(., depvar) X_raw = st_data(., tokens(indepvars)) panel = st_data(., panelid) K = cols(X_raw) // Reshape Y to T x N (each column = one unit's time series) Y_wide = J(T, N, .) for (i = 1; i <= N; i++) { Y_wide[., i] = Y_raw[(i-1)*T+1 :: i*T] } // GAUSS cum() lines 13-15 // t = rows(y); n = cols(y); k = cols(x)/n // Here: t = T, n = N, k = K // GAUSS lines 18-22: prepend zeros then difference // de = zeros(t, n) // y = zeros(1,n)|y => (T+1) x N // x = zeros(1,k*n)|x => (T+1) x (k*N) // dy = diff(y,1) => T x N // dx = diff(x,1) => T x (k*N) de = J(T, N, 0) for (i = 1; i <= N; i++) { // Extract unit i's data real colvector yi_level, dyi real matrix xi_level, dxi real scalar r_start, r_end yi_level = Y_wide[., i] r_start = (i-1)*T + 1 r_end = i*T xi_level = X_raw[r_start::r_end, .] // GAUSS zero-padding + differencing (lines 19-22): // Prepend zero row, then take first differences // Result: dy[1] = y[1] - 0 = y[1], dy[2] = y[2]-y[1], ... // This gives T observations dyi = yi_level - (J(1, 1, 0) \ yi_level[1::T-1]) dxi = xi_level - (J(1, K, 0) \ xi_level[1::T-1, .]) // GAUSS line 27: de[.,i] = (eye(t) - dxi*inv(dxi'dxi)*dxi')*dyi if (predet != 0) { // Predetermined coefficient: de = dy - predet*dx de[., i] = dyi - dxi * J(K, 1, predet) } else { // OLS projection: de = (I - X(X'X)^{-1}X') * dy de[., i] = dyi - dxi * (invsym(cross(dxi, dxi)) * cross(dxi, dyi)) } } // GAUSS line 32: nf = fact(de, kmax, pen, cri) nf = _xtdh_fact(de, kmax, pen, cri) // GAUSS line 33: {f, lam} = prin(de, nf) _xtdh_prin(de, nf, F, Lambda) // GAUSS line 34: de = de - f*lam' de = de - F * Lambda' // GAUSS line 35: e = cumsumc(de) // Cumulative sum column-by-column e = J(T, N, 0) for (i = 1; i <= N; i++) { e[., i] = runningsum(de[., i]) } } /* ** Individual DH statistic ** Replicates GAUSS gdh() lines 40-63 EXACTLY ** ** GAUSS code: ** wl = w[1:t-1] ** w0 = w[2:t] ** b1 = inv(wl'wl)*(wl'w0) // OLS phi_hat ** b2 = inv(wl'w0)*(w0'w0) // IV phi_tilde ** e = w0 - wl*b1 ** s = (e'e)/rows(e) // sigma^2 ** io = fejer(e,p) // one-sided LR cov ** sig = s + io + io' // omega^2 ** s = s^2/sig // sigma^4/omega^2 ** dhs = (b2-b1)^2/(s*inv(wl'wl)) // = (omega^2/sigma^4)*(b2-b1)^2*(wl'wl) */ real scalar _xtdh_gdh(real colvector w, real scalar p) { real scalar t, b1, b2, s, io, sig, dhs real colvector wl, w0, resid t = rows(w) // GAUSS lines 45-46 wl = w[1::t-1] w0 = w[2::t] // GAUSS line 48: b1 = inv(wl'wl)*(wl'w0) [OLS] b1 = (wl' * w0) / (wl' * wl) // GAUSS line 49: b2 = inv(wl'w0)*(w0'w0) [IV] b2 = (w0' * w0) / (wl' * w0) // GAUSS line 50: e = w0 - wl*b1 resid = w0 - wl * b1 // GAUSS line 51: s = (e'e)/rows(e) => sigma^2 s = (resid' * resid) / rows(resid) // GAUSS lines 52-56: one-sided long-run covariance if (p == 0) { io = 0 } else { io = _xtdh_fejer(resid, p) } // GAUSS line 57: sig = s + io + io' => omega^2 = sigma^2 + 2*io // (io is scalar, io' = io) sig = s + 2 * io // GAUSS line 58: s = s^2/sig => sigma^4 / omega^2 s = s^2 / sig // GAUSS line 60: dhs = (b2-b1)^2/(s*inv(wl'wl)) // = (b2-b1)^2 / ((sigma^4/omega^2) * (1/(wl'wl))) // = (b2-b1)^2 * (wl'wl) * (omega^2/sigma^4) // = Si * (phi_tilde - phi_hat)^2 * sum(e_{t-1}^2) dhs = (b2 - b1)^2 / (s * (1 / (wl' * wl))) return(dhs) } /* ** Panel DH statistic ** Replicates GAUSS pdh() lines 65-102 EXACTLY ** ** Note: GAUSS line 71 drops first column: w = w[., 2:cols(w)] ** This is application-specific for the Fisher effect analysis. ** For a general-purpose Stata command, we keep all columns. */ real scalar _xtdh_pdh(real matrix w, real scalar p) { real scalar t, n, i, b1, b2, s, io, snt, lnt, dhs real colvector wl_pool, w0_pool, resid_i real matrix wl, w0, e_reshaped t = rows(w) n = cols(w) // GAUSS lines 72-73: create lags wl = w[1::t-1, .] w0 = w[2::t, .] // GAUSS lines 74-75: pool via vec (column stacking) wl_pool = vec(wl) w0_pool = vec(w0) // GAUSS line 77: b1 = inv(wl'wl)*(wl'w0) [pooled OLS] b1 = (wl_pool' * w0_pool) / (wl_pool' * wl_pool) // GAUSS line 78: b2 = inv(wl'w0)*(w0'w0) [pooled IV] b2 = (w0_pool' * w0_pool) / (wl_pool' * w0_pool) // GAUSS line 83: e = reshape((w0-wl*b1), n, t-1)' // GAUSS reshape(vec, nrows, ncols) fills row-by-row // vec(wl) stacks columns, so (w0_pool - wl_pool*b1) is n*(t-1) x 1 // reshape to n x (t-1), then transpose to (t-1) x n e_reshaped = colshape(w0_pool - wl_pool * b1, n) // GAUSS lines 80-97: compute variance components snt = 0 // sum of sigma_i^2 lnt = 0 // sum of omega_i^2 for (i = 1; i <= n; i++) { resid_i = e_reshaped[., i] // GAUSS line 86: s = (e[.,i]'e[.,i])/rows(e) s = (resid_i' * resid_i) / rows(resid_i) // GAUSS lines 87-91: one-sided long-run covariance if (p == 0) { io = 0 } else { io = _xtdh_fejer(resid_i, p) } // GAUSS lines 93-94 snt = snt + s // sum sigma_i^2 lnt = lnt + (s + 2 * io) // sum omega_i^2 (io+io' = 2*io for scalar) } // GAUSS line 99: // dhs = ((lnt/n)*(b2-b1)^2) / (((snt/n)^2)*inv(wl'wl)) // = Sn * (b2-b1)^2 * (wl'wl) // where Sn = (lnt/n) / (snt/n)^2 = omega_bar^2 / sigma_bar^4 dhs = ((lnt / n) * (b2 - b1)^2) / (((snt / n)^2) * (1 / (wl_pool' * wl_pool))) return(dhs) } /* ** Factor number selection ** Replicates GAUSS fact() lines 104-158 EXACTLY */ real scalar _xtdh_fact(real matrix e, real scalar kmax, real scalar p, real scalar c) { real scalar T_f, N_f, k, s, smax, pen_val real colvector cr real matrix F, Lambda, u T_f = rows(e) N_f = cols(e) // GAUSS lines 109-111: compute maximum variance _xtdh_prin(e, kmax, F, Lambda) u = e - F * Lambda' smax = sum(u:^2) / (N_f * T_f) cr = J(kmax, 1, .) // GAUSS lines 113-153 for (k = 1; k <= kmax; k++) { _xtdh_prin(e, k, F, Lambda) u = e - F * Lambda' s = sum(u:^2) / (N_f * T_f) // GAUSS lines 121-140: penalty selection if (c == 1 | c == 2) { if (p == 1) { pen_val = (N_f + T_f) / (N_f * T_f) * ln((N_f * T_f) / (N_f + T_f)) } else if (p == 2) { pen_val = (N_f + T_f) / (N_f * T_f) * ln(min((N_f, T_f))) } else { pen_val = ln(min((N_f, T_f))) / min((N_f, T_f)) } } else { // c == 3 (AIC/BIC) if (p == 1) { pen_val = 2 * (N_f + T_f - k) / (N_f * T_f) } else { pen_val = (N_f + T_f - k) * ln(N_f * T_f) / (N_f * T_f) } } // GAUSS lines 142-150: criterion value if (c == 1) { cr[k] = s + k * smax * pen_val // PC } else if (c == 2) { cr[k] = ln(s) + k * pen_val // IC } else { cr[k] = s + k * smax * pen_val // AIC/BIC } } // GAUSS lines 155-157: return k that minimizes criterion // GAUSS: cr = sortc(seqa(1,1,rows(cr))~cr, 2); retp(cr[1,1]) return(selectindex(cr :== min(cr))[1]) } /* ** Principal components estimation ** Replicates GAUSS prin() lines 160-178 EXACTLY ** ** GAUSS code: ** if n > t: ** {f0, v, f} = svd1(e*e') ** f = f0[.,1:nf]*sqrt(t) ** lam = (e'f)/t ** else: ** {f0, v, f} = svd1(e'e) ** lam = f0[.,1:nf]*sqrt(n) ** f = (e*lam)/n */ void _xtdh_prin(real matrix e, real scalar nf, real matrix F, real matrix Lambda) { real scalar T_p, N_p real matrix U, Vt real colvector sv T_p = rows(e) N_p = cols(e) // GAUSS lines 166-175 if (N_p > T_p) { // svd(e*e') => eigendecomposition of T x T matrix svd(e * e', U, sv, Vt) F = U[., 1::nf] * sqrt(T_p) Lambda = (e' * F) / T_p } else { // svd(e'e) => eigendecomposition of N x N matrix svd(e' * e, U, sv, Vt) Lambda = U[., 1::nf] * sqrt(N_p) F = (e * Lambda) / N_p } } /* ** Fejer (Bartlett) kernel for one-sided long-run covariance ** ** Computes: io = sum_{j=1}^{M} w_j * gamma_j ** where: ** w_j = 1 - j/(M+1) (Bartlett kernel weight) ** gamma_j = (1/T) * sum_{t=j+1}^{T} e_t * e_{t-j} (autocovariance) ** ** This is the ONE-SIDED covariance. ** Full LR variance: omega^2 = sigma^2 + 2*io */ real scalar _xtdh_fejer(real colvector e, real scalar M) { real scalar T_k, j, io, w, gam_j T_k = rows(e) io = 0 for (j = 1; j <= M; j++) { w = 1 - j / (M + 1) gam_j = (e[1::T_k-j]' * e[j+1::T_k]) / T_k io = io + w * gam_j } return(io) } end