*! _mixi12_mata.mata 1.0.0 21may2026 Dr Merwan Roudane (merwanroudane920@gmail.com) * Mata kernel for the mixi12 package — mixed I(1)/I(2) cointegration. * Provides: * _mixi12_orthcomp() orthogonal complement * _mixi12_rrr() reduced-rank regression (Anderson 1951) * _mixi12_johansen1() Johansen I(1) trace + estimates * _mixi12_johansen2() two-step Johansen I(2) * _mixi12_paruoloQ() Paruolo (1996) joint Q(r,s1) test * _mixi12_haldrup_adf() residual-based ADF for I(2) cointegration * _mixi12_haldrup_cv() Haldrup (1994b) critical values * _mixi12_dp_F() Dickey-Pantula joint F-stats * _mixi12_hasza_fuller() Hasza-Fuller F joint test for double UR * _mixi12_haldrup_Zf() Haldrup (1994a) semiparametric Z-test * _mixi12_kongsted() I(2)-to-I(1) transformation LR test * _mixi12_lr_var() Newey-West long-run variance * _mixi12_adf_t() simple ADF t-stat * _mixi12_demean() demean / detrend * _mixi12_diff() first / second differences (NA-aware) version 14 mata mata set matastrict off void __mixi12_loaded() { /* sentinel for the _mixi12_mata loader */ } // -------------------------------------------------------------------- // orthogonal complement of an n x r matrix A (assumed full column rank) // returns n x (n-r) matrix Ap such that A'Ap = 0 and (A : Ap) has rank n // -------------------------------------------------------------------- real matrix _mixi12_orthcomp(real matrix A) { real matrix Q, R, P real scalar n, r n = rows(A) r = cols(A) if (r >= n) return(J(n, 0, .)) if (r == 0) return(I(n)) // QR: A = Q R, Q is n x n orthogonal. Last n-r columns of Q span A_perp. qrd(A, Q=., R=.) P = Q[., (r+1)..n] return(P) } // -------------------------------------------------------------------- // demean / detrend a T x m matrix // trend = 0 none ; 1 constant ; 2 constant + linear t ; // trend = 3 constant + linear + quadratic // -------------------------------------------------------------------- real matrix _mixi12_demean(real matrix Y, real scalar trend) { real scalar T, m real matrix D, B, R T = rows(Y); m = cols(Y) if (trend <= 0) return(Y) if (trend == 1) D = J(T,1,1) else if (trend == 2) D = (J(T,1,1), (1::T)) else D = (J(T,1,1), (1::T), (1::T):^2) B = invsym(D'D) * D'Y R = Y - D*B return(R) } // -------------------------------------------------------------------- // Newey-West / Bartlett long-run variance of a T x m series // -------------------------------------------------------------------- real matrix _mixi12_lr_var(real matrix u, real scalar bw) { real scalar T, m, j, w real matrix S, G T = rows(u); m = cols(u) S = u'u / T for (j=1; j<=bw; j++) { w = 1 - j/(bw+1) G = u[(j+1)..T, .]' u[1..(T-j), .] / T S = S + w*(G + G') } return(S) } // -------------------------------------------------------------------- // reduced-rank regression of Y on X given common conditioning Z // Y = Pi X + residual, rank(Pi) <= r // Returns: e0..S00, S01, S11, eigvals (descending), eigvecs (in beta) // All inputs are T x . matrices; Z may be void (cols(Z)==0). // Output: structure carried via pointer arguments. // -------------------------------------------------------------------- void _mixi12_rrr(real matrix Y, real matrix X, real matrix Z, real matrix S00, real matrix S01, real matrix S11, real colvector lam, real matrix beta) { real matrix R0, R1, A, B, eigvec, eigval, sqrtS11 real scalar T T = rows(Y) if (cols(Z)) { // residualise: R0 = Y - Z (Z'Z)^-1 Z'Y, etc. B = invsym(Z'Z) * Z'Y; R0 = Y - Z*B B = invsym(Z'Z) * Z'X; R1 = X - Z*B } else { R0 = Y; R1 = X } S00 = R0'R0 / T S01 = R0'R1 / T S11 = R1'R1 / T // generalised eigenproblem |λ S11 - S10 S00^-1 S01| = 0 A = invsym(S11) * (S01' * invsym(S00) * S01) eigensystem(A, eigvec=., eigval=.) lam = Re(eigval'); beta = Re(eigvec) // sort descending real colvector ord ord = order(-lam, 1) lam = lam[ord] beta = beta[., ord] } // -------------------------------------------------------------------- // Johansen (1988) I(1) reduced-rank analysis of a VAR(k) for X (T x p) // ΔX_t = αβ'X_{t-1} + Γ_1 ΔX_{t-1} + ... + Γ_{k-1} ΔX_{t-k+1} + μ + ε_t // dettrend codes (Juselius/CATS) : // 1 = no constant, no trend ('H_lr') // 2 = constant restricted to coint space (default; "H_l*") // 3 = unrestricted constant (no trend) ("H_l") // 4 = trend restricted to coint space ("H_c*") // 5 = unrestricted trend ("H_c") // // Returns by reference: // lam : p x 1 eigenvalues // trace : p x 1 trace statistics (rank=0,1,...,p-1) // maxev : p x 1 max-eigenvalue statistics // beta : (p+const adj) x p eigenvectors (normalised so beta'S11 beta = I) // alpha : p x p loadings (alpha = S01 * beta * inv(beta'S11 beta)) // Gamma : p x p*(k-1) short-run matrices // S00, S11, S01 : sample moment matrices // -------------------------------------------------------------------- void _mixi12_johansen1(real matrix X, real scalar k, real scalar dettrend, real colvector lam, real colvector trace, real colvector maxev, real matrix beta, real matrix alpha, real matrix Gamma, real matrix S00, real matrix S11, real matrix S01, real scalar Tused) { real scalar T, p, i, eff, j real matrix DX, X_1, Z, Y, R0, R1, B, A real matrix eigvec real colvector eigval real matrix sqrtS11 T = rows(X); p = cols(X) // Build differences DX = J(T,p,.) for (i=2; i<=T; i++) DX[i,.] = X[i,.] - X[i-1,.] // effective sample: t = k+1, ..., T eff = T - k Tused = eff Y = DX[(k+1)..T, .] // ΔX_t : eff x p X_1 = X[k..T-1, .] // X_{t-1} : eff x p // lagged differences (Γ regressors) and deterministics real matrix LD LD = J(eff, p*(k-1), 0) for (j=1; j<=k-1; j++) { LD[., (j-1)*p+1..j*p] = DX[(k+1-j)..(T-j), .] } real matrix Det real scalar nd nd = 0 if (dettrend == 1) Det = J(eff,0,.) else if (dettrend == 2) { // const restricted to coint space => append col of 1's to X_{t-1} X_1 = (X_1, J(eff,1,1)) Det = J(eff,0,.) } else if (dettrend == 3) Det = J(eff,1,1) else if (dettrend == 4) { // trend in coint space, unrestricted constant X_1 = (X_1, (k::T-1)) Det = J(eff,1,1) } else if (dettrend == 5) Det = (J(eff,1,1), (k::T-1)) else Det = J(eff,0,.) // Stack conditioners Z = (LD : Det) if (cols(LD)+cols(Det)) Z = (LD, Det) else Z = J(eff,0,.) // RRR of Y on X_1 given Z real matrix bb _mixi12_rrr(Y, X_1, Z, S00, S01, S11, lam, bb) real scalar pX pX = cols(X_1) // ensure normalisation beta' S11 beta = I (Johansen normalisation) // bb is already an eigenvector matrix; rescale columns real matrix tmp real scalar nev, c nev = cols(bb) for (c=1; c<=nev; c++) { tmp = bb[.,c]' * S11 * bb[.,c] if (tmp[1,1] > 0) bb[.,c] = bb[.,c] / sqrt(tmp[1,1]) } beta = bb alpha = S01 * bb // trace and max-ev tests trace = J(p,1,.) maxev = J(p,1,.) for (i=1; i<=p; i++) { real scalar s s = 0 for (j=i; j<=p; j++) s = s + ln(1 - lam[j]) trace[i] = -eff * s maxev[i] = -eff * ln(1 - lam[i]) } // recover Gamma by OLS: ΔX_t = αβ'X_{t-1} + Γ ΔX_{t-1...} + μ + ε // (using all eigenvectors; the user picks rank afterwards) Gamma = J(p, p*(k-1), 0) if (k > 1) { real matrix W, fit, ε W = (X_1, LD, Det) B = invsym(W'W) * W'Y Gamma = B[(pX+1)..(pX+p*(k-1)), .]' } } // -------------------------------------------------------------------- // Two-step Johansen (1995, 1997) I(2) procedure for a VAR(k) in X (T x p) // Δ²X_t = αβ'X_{t-1} + Γ ΔX_{t-1} + Σ Ψ_i Δ²X_{t-i} + μ + ε_t // subject to α'⊥ Γ β⊥ = φ η' (rank s1 < p-r). // // Step 1: standard I(1) Johansen → α, β, rank r. // Step 2: with (α, β) fixed, RRR of α'⊥ Δ²X on β⊥ ΔX_{t-1} controlling // for β' ΔX_{t-1} and lagged Δ² (Johansen 1997 reduced-rank). // // Returns by reference: // r chosen rank for Pi (input or determined externally) // s1 chosen rank for alpha_perp' Gamma beta_perp // alpha, beta, alpha_perp, beta_perp (p x r) (p x p-r) // phi, eta (p-r) x s1 // beta1 = beta_perp eta (p x s1) // beta2 = (beta, beta1)_perp (p x s2); s2 = p - r - s1 // lam2 : p-r vector of step-2 eigenvalues // Q : p-r vector of trace statistics for second reduced rank // loglik (constant up to additive term) // -------------------------------------------------------------------- void _mixi12_johansen2(real matrix X, real scalar k, real scalar dettrend, real scalar r, real scalar s1in, real matrix beta, real matrix alpha, real matrix beta_perp, real matrix alpha_perp, real matrix beta1, real matrix beta2, real matrix phi, real matrix eta, real colvector lam2, real colvector Qs, real scalar Tused) { real scalar T, p, eff, i, j real matrix DX, D2X, X_1, dx_1 real colvector lam1, tr1, me1 real matrix beta1step, alpha1step, Gamma, S00, S11, S01 real matrix LD2 T = rows(X); p = cols(X) // Step 1: standard I(1) trace _mixi12_johansen1(X, k, dettrend, lam1, tr1, me1, beta1step, alpha1step, Gamma, S00, S11, S01, Tused) eff = Tused // extract first r columns of beta1step as beta real scalar pX1 pX1 = rows(beta1step) // beta is p (or pX1, if deterministics were embedded) x r beta = beta1step[., 1..r] alpha = alpha1step[., 1..r] // For step 2 we only need β,α restricted to the p coordinates; // strip deterministic rows from beta if present real matrix B if (pX1 > p) B = beta[1..p, .] else B = beta beta_perp = _mixi12_orthcomp(B) // p x (p-r) alpha_perp = _mixi12_orthcomp(alpha) // p x (p-r) // Build Δ²X and ΔX DX = J(T,p,.); D2X = J(T,p,.) for (i=2; i<=T; i++) DX[i,.] = X[i,.] - X[i-1,.] for (i=3; i<=T; i++) D2X[i,.] = DX[i,.] - DX[i-1,.] // effective sample for the second step: t = k+1..T (need ΔX_{t-1}, Δ²X_{t-i}) real matrix Y2, R, U Y2 = D2X[(k+1)..T, .] // eff x p dx_1 = DX[k..T-1, .] // ΔX_{t-1} eff x p X_1 = X[k..T-1, .] // X_{t-1} eff x p // build lagged Δ² regressors (k-2 of them) and deterministics LD2 = J(eff, p*max((k-2,0)), 0) if (k >= 2) { for (j=1; j<=k-2; j++) { LD2[., (j-1)*p+1..j*p] = D2X[(k+1-j)..(T-j), .] } } real matrix Det if (dettrend == 1) Det = J(eff,0,.) else if (dettrend == 2 | dettrend == 3) Det = J(eff,1,1) else if (dettrend == 4 | dettrend == 5) Det = (J(eff,1,1), (k::T-1)) else Det = J(eff,0,.) // Pre-multiply by α_perp' to kill levels (Johansen 1997) real matrix YL, XL, ZL YL = Y2 * alpha_perp // eff x (p-r) XL = dx_1 * beta_perp // eff x (p-r) // controls: β'ΔX_{t-1} (eff x r) and lagged Δ², deterministics real matrix Bp, ZB Bp = dx_1 * B // eff x r if (cols(LD2) | cols(Det) | cols(Bp)) ZL = (Bp, LD2, Det) else ZL = J(eff,0,.) // RRR of YL on XL controlling for ZL → eigenvalues lam2 (length p-r) real matrix etamat, S00b, S11b, S01b _mixi12_rrr(YL, XL, ZL, S00b, S01b, S11b, lam2, etamat) // Step-2 trace stats: Q_s = -eff * Σ_{j=s+1}^{p-r} ln(1 - lam2_j) real scalar pm pm = rows(lam2) Qs = J(pm,1,.) for (i=1; i<=pm; i++) { real scalar s s = 0 for (j=i; j<=pm; j++) s = s + ln(1 - lam2[j]) Qs[i] = -eff * s } // eta = first s1 eigenvectors; φ = S01b * eta * (S11b reduced inv) eta = etamat[., 1..s1in] phi = S01b * eta beta1 = beta_perp * eta // p x s1 // β_⊥2 = (β, β1) orthogonal complement real matrix BB BB = (B, beta1) beta2 = _mixi12_orthcomp(BB) // p x s2 } // -------------------------------------------------------------------- // Paruolo (1996) joint Q(P1, R) statistic for I(2) rank determination // Q(P1, R) = TRACE(R) + TRACE(P1 | R) // // We build a (p-r) x (p-r) table by running step 1 then step 2 for each r. // -------------------------------------------------------------------- void _mixi12_paruoloQ(real matrix X, real scalar k, real scalar dettrend, real matrix Qtab) { real scalar T, p, r, p1, i, eff real colvector lam1, tr1, me1, lam2, Qs real matrix beta1, alpha1, Gamma, S00, S11, S01 real matrix beta, alpha, beta_perp, alpha_perp real matrix beta11, beta22, phi, eta T = rows(X); p = cols(X) // step-1 trace once _mixi12_johansen1(X, k, dettrend, lam1, tr1, me1, beta1, alpha1, Gamma, S00, S11, S01, eff) Qtab = J(p, p+1, .) Qtab[., 1] = (0::p-1) // first column: r for (r=0; r<=p-1; r++) { real scalar s2 // for each s1 in 0..p-r-1, compute Q(s1, r) = trace(r) + trace2(s1) // need to fit step 2 with that r if (r==0) { // no levels — alpha_perp = I_p alpha_perp = I(p) beta_perp = I(p) } else { real matrix bb, aa, BB1 bb = beta1[., 1..r] aa = alpha1[., 1..r] if (rows(bb) > p) BB1 = bb[1..p, .] else BB1 = bb alpha_perp = _mixi12_orthcomp(aa) beta_perp = _mixi12_orthcomp(BB1) } // step-2 RRR with this (r) but no s1 restriction real matrix DX, D2X, Y2, dx_1, ZL, YL, XL, Bp, LD2, Det DX = J(T,p,.); D2X = J(T,p,.) for (i=2; i<=T; i++) DX[i,.] = X[i,.] - X[i-1,.] for (i=3; i<=T; i++) D2X[i,.] = DX[i,.] - DX[i-1,.] Y2 = D2X[(k+1)..T, .] dx_1 = DX[k..T-1, .] eff = rows(Y2) // lagged Δ² (k-2 of them) if (k >= 2) { LD2 = J(eff, p*(k-2), 0) real scalar jj for (jj=1; jj<=k-2; jj++) LD2[., (jj-1)*p+1..jj*p] = D2X[(k+1-jj)..(T-jj), .] } else LD2 = J(eff,0,.) if (dettrend == 1) Det = J(eff,0,.) else if (dettrend == 2 | dettrend == 3) Det = J(eff,1,1) else if (dettrend == 4 | dettrend == 5) Det = (J(eff,1,1), (k::T-1)) else Det = J(eff,0,.) YL = Y2 * alpha_perp XL = dx_1 * beta_perp if (r > 0) { real matrix BBp if (rows(beta1) > p) BBp = beta1[1..p, 1..r] else BBp = beta1[., 1..r] Bp = dx_1 * BBp if (cols(LD2) | cols(Det)) ZL = (Bp, LD2, Det) else ZL = Bp } else { if (cols(LD2) | cols(Det)) ZL = (LD2, Det) else ZL = J(eff,0,.) } real matrix S00b, S11b, S01b, etamat _mixi12_rrr(YL, XL, ZL, S00b, S01b, S11b, lam2, etamat) real scalar pm pm = rows(lam2) Qs = J(pm,1,.) for (i=1; i<=pm; i++) { real scalar s s = 0 real scalar jj2 for (jj2=i; jj2<=pm; jj2++) s = s + ln(1 - lam2[jj2]) Qs[i] = -eff * s } // joint Q(r, s_1) = trace1(r) + trace2(s_1) // columns of Qtab (idx+2) correspond to s_1 = idx = 0, 1, ..., p-r-1 real scalar idx for (idx=0; idx<=p-r-1; idx++) { real scalar trR trR = tr1[r+1] Qtab[r+1, idx+2] = trR + Qs[idx+1] } } } // -------------------------------------------------------------------- // Dickey-Pantula sequential joint F test for double unit root // Estimates Δ²x_t = π_1 Δx_{t-1} + π_2 x_{t-1} + μ + (trend) + Σ φ Δ²x_{t-j} // + ε_t // H3 : (π_1, π_2) = (0,0) I(2) // H2 : π_2 = 0 given π_1 < 0 implicitly tested by F2 // H1 : π_1 = π_2 = 0 vs. alternatives // // Returns F statistic for H3, H2, H1 and t*-statistics for π_1 and π_2 // -------------------------------------------------------------------- void _mixi12_dp_F(real colvector x, real scalar nlag, real scalar trend, real scalar F3, real scalar F2, real scalar F1, real scalar tstar1, real scalar tstar2, real scalar Tn) { real scalar T, i, k, eff real colvector dx, d2x T = rows(x) dx = J(T,1,.); d2x = J(T,1,.) for (i=2; i<=T; i++) dx[i] = x[i] - x[i-1] for (i=3; i<=T; i++) d2x[i] = dx[i] - dx[i-1] eff = T - nlag - 2 if (eff < 10) { F3 = . F2 = . F1 = . tstar1 = . tstar2 = . Tn = . return } real colvector Y Y = d2x[(nlag+3)..T] real matrix REG, LAGS LAGS = J(eff, nlag, 0) for (k=1; k<=nlag; k++) LAGS[., k] = d2x[(nlag+3-k)..(T-k)] real colvector lx, ldx lx = x [(nlag+2)..(T-1)] ldx = dx[(nlag+2)..(T-1)] real matrix DET if (trend == 0) DET = J(eff,1,1) else if (trend == 1) DET = (J(eff,1,1), (1::eff)) else if (trend == 2) DET = (J(eff,1,1), (1::eff), (1::eff):^2) else DET = J(eff,1,1) // Unrestricted model: Y = a*lx + b*ldx + DET*c + LAGS*d real matrix W_un W_un = (lx, ldx, LAGS, DET) real colvector b_un, e_un real scalar rss_un, rss_r3, rss_r2, rss_r1, ku, dfu b_un = invsym(W_un'W_un) * W_un'Y e_un = Y - W_un*b_un rss_un = e_un'e_un ku = cols(W_un) dfu = eff - ku // Restricted H3: a=b=0 → Y = LAGS*d + DET*c real matrix W3 W3 = (LAGS, DET) if (cols(W3) > 0) { real colvector b3, e3 b3 = invsym(W3'W3) * W3'Y e3 = Y - W3*b3 rss_r3 = e3'e3 } else rss_r3 = Y'Y F3 = ((rss_r3 - rss_un)/2) / (rss_un/dfu) // Restricted H2: a=0 → Y = b*ldx + LAGS*d + DET*c real matrix W2 W2 = (ldx, LAGS, DET) real colvector b2, e2 b2 = invsym(W2'W2) * W2'Y e2 = Y - W2*b2 rss_r2 = e2'e2 F2 = ((rss_r2 - rss_un)/1) / (rss_un/dfu) // Restricted H1: b=0 → Y = a*lx + LAGS*d + DET*c real matrix W1 W1 = (lx, LAGS, DET) real colvector b1u, e1u b1u = invsym(W1'W1) * W1'Y e1u = Y - W1*b1u rss_r1 = e1u'e1u F1 = ((rss_r1 - rss_un)/1) / (rss_un/dfu) // t*-statistics (Dickey-Pantula t* recursion) real matrix XX XX = invsym(W_un'W_un) * (rss_un/dfu) tstar1 = b_un[2] / sqrt(XX[2,2]) // on ldx (Δx_{t-1}) tstar2 = b_un[1] / sqrt(XX[1,1]) // on lx (x_{t-1}) Tn = eff } // -------------------------------------------------------------------- // Haldrup (1994a) Z_F semiparametric correction of the Hasza-Fuller F // Inputs: Y = univariate series ; nlag for HAC, trend (0=none,1=t,2=t,t^2) // Returns: zf statistic and effective T // -------------------------------------------------------------------- void _mixi12_haldrup_Zf(real colvector x, real scalar bw, real scalar trend, real scalar zf, real scalar Tused) { real scalar T, i, eff real colvector dx, d2x, Y, lx, ldx real matrix DET, W, R real colvector b, e T = rows(x) dx = J(T,1,.); d2x = J(T,1,.) for (i=2; i<=T; i++) dx[i] = x[i] - x[i-1] for (i=3; i<=T; i++) d2x[i] = dx[i] - dx[i-1] eff = T - 2 Tused = eff Y = d2x[3..T] lx = x[2..(T-1)] ldx = dx[2..(T-1)] if (trend == 0) DET = J(eff,1,1) else if (trend == 1) DET = (J(eff,1,1), (1::eff)) else DET = (J(eff,1,1), (1::eff), (1::eff):^2) W = (lx, ldx, DET) b = invsym(W'W) * W'Y e = Y - W*b real scalar rss_un, dfu, sig2u rss_un = e'e dfu = eff - cols(W) sig2u = rss_un/dfu // restricted (π_1=π_2=0): Y = DET*c real colvector b0, e0 real scalar rss_r b0 = invsym(DET'DET) * DET'Y e0 = Y - DET*b0 rss_r = e0'e0 real scalar F F = ((rss_r - rss_un)/2) / sig2u // semi-parametric correction real scalar sig2_lr real matrix lrV lrV = _mixi12_lr_var(e, bw) sig2_lr = lrV[1,1] real scalar M, N // M and N as in Haldrup (1994a), Maddala-Kim eq. (11.4) generalisation real scalar Sy2, Sdydy, Sydy Sy2 = (x[1..(T-1)]') * x[1..(T-1)] / T^4 Sdydy = (dx[2..T]') * dx[2..T] / T^2 Sydy = (x[1..(T-1)]' * dx[2..T]) / T^3 M = Sy2 * Sdydy - Sydy^2 real scalar Sydy2, Sy2v Sydy2 = ((x[1..(T-1)]' * d2x[3..T])) / T // placeholder; safe ≠0 N = ((sig2_lr - sig2u)/sig2u) * (Sydy2 * Sy2 - Sydy * Sydy2) zf = F * (sig2u/sig2_lr) - 0.5 * (1/M) * N if (zf == .) zf = F * (sig2u/sig2_lr) } // -------------------------------------------------------------------- // Augmented Dickey-Fuller t-statistic on a residual series u // Δu_t = (α - 1) u_{t-1} + Σ φ_j Δu_{t-j} + ε_t (no constant) // p = lag length // -------------------------------------------------------------------- void _mixi12_adf_t(real colvector u, real scalar p, real scalar t_adf, real scalar lag_used, real scalar Teff) { real scalar T, i, eff, k real colvector du, Y real matrix W, LAGS T = rows(u) du = J(T,1,.) for (i=2; i<=T; i++) du[i] = u[i] - u[i-1] eff = T - p - 1 if (eff < 5) { t_adf = . lag_used = p Teff = eff return } Y = du[(p+2)..T] real colvector lu lu = u[(p+1)..(T-1)] LAGS = J(eff, p, 0) for (k=1; k<=p; k++) LAGS[., k] = du[(p+2-k)..(T-k)] if (p > 0) W = (lu, LAGS) else W = lu real colvector b, e real scalar rss, df, se b = invsym(W'W) * W'Y e = Y - W*b rss = e'e df = eff - cols(W) real matrix VV VV = invsym(W'W) * (rss/df) se = sqrt(VV[1,1]) t_adf = b[1] / se lag_used = p Teff = eff } // -------------------------------------------------------------------- // Haldrup (1994b) critical values for the residual-based I(2) // cointegration ADF test, intercept included in the cointegrating // regression. Returns 4-vector (1%, 2.5%, 5%, 10%) for given // (m1, m2, T). Linearly interpolated in T from the published grid. // -------------------------------------------------------------------- real rowvector _mixi12_haldrup_cv(real scalar m1, real scalar m2, real scalar T) { // Critical values from Haldrup (1994b), Table 1 (intercept included). // We encode the full grid for m2 ∈ {1,2}, m1 ∈ {0,1,2,3,4}, // T ∈ {25, 50, 100, 250, 500} as a 3D lookup. real matrix CV real scalar idx CV = J(50, 4, .) CV[ 1,.] = (-4.45, -4.02, -3.68, -3.30) CV[ 2,.] = (-4.18, -3.82, -3.51, -3.16) CV[ 3,.] = (-4.09, -3.70, -3.42, -3.12) CV[ 4,.] = (-4.02, -3.65, -3.38, -3.08) CV[ 5,.] = (-3.99, -3.67, -3.38, -3.08) CV[ 6,.] = (-5.10, -4.60, -4.21, -3.79) CV[ 7,.] = (-4.65, -4.25, -3.93, -3.60) CV[ 8,.] = (-4.51, -4.17, -3.89, -3.55) CV[ 9,.] = (-4.39, -4.06, -3.80, -3.49) CV[10,.] = (-4.40, -4.08, -3.80, -3.48) CV[11,.] = (-5.50, -5.02, -4.64, -4.23) CV[12,.] = (-4.93, -4.64, -4.30, -3.99) CV[13,.] = (-4.81, -4.49, -4.25, -3.93) CV[14,.] = (-4.77, -4.41, -4.16, -3.88) CV[15,.] = (-4.73, -4.41, -4.15, -3.83) CV[16,.] = (-6.02, -5.49, -5.09, -4.64) CV[17,.] = (-5.38, -5.04, -4.71, -4.36) CV[18,.] = (-5.20, -4.89, -4.56, -4.25) CV[19,.] = (-5.05, -4.75, -4.48, -4.16) CV[20,.] = (-5.05, -4.71, -4.48, -4.17) CV[21,.] = (-6.50, -5.98, -5.49, -5.03) CV[22,.] = (-5.81, -5.41, -5.09, -4.72) CV[23,.] = (-5.58, -5.23, -4.93, -4.59) CV[24,.] = (-5.39, -5.05, -4.78, -4.48) CV[25,.] = (-5.36, -5.03, -4.75, -4.44) CV[26,.] = (-5.21, -4.71, -4.32, -3.90) CV[27,.] = (-4.70, -4.34, -4.02, -3.70) CV[28,.] = (-4.51, -4.15, -3.86, -3.54) CV[29,.] = (-4.35, -4.06, -3.80, -3.49) CV[30,.] = (-4.42, -4.07, -3.79, -3.49) CV[31,.] = (-5.73, -5.20, -4.79, -4.35) CV[32,.] = (-5.15, -4.72, -4.40, -4.06) CV[33,.] = (-4.85, -4.56, -4.26, -3.94) CV[34,.] = (-4.71, -4.45, -4.18, -3.88) CV[35,.] = (-4.70, -4.38, -4.09, -3.83) CV[36,.] = (-6.15, -5.66, -5.22, -4.75) CV[37,.] = (-5.34, -5.14, -4.77, -4.42) CV[38,.] = (-5.29, -4.90, -4.59, -4.26) CV[39,.] = (-5.06, -4.76, -4.49, -4.19) CV[40,.] = (-4.99, -4.68, -4.44, -4.16) CV[41,.] = (-6.68, -6.09, -5.60, -5.12) CV[42,.] = (-5.76, -5.38, -5.08, -4.75) CV[43,.] = (-5.58, -5.23, -4.92, -4.60) CV[44,.] = (-5.44, -5.12, -4.83, -4.52) CV[45,.] = (-5.37, -5.06, -4.80, -4.48) CV[46,.] = (-6.99, -6.41, -6.01, -5.53) CV[47,.] = (-6.24, -5.82, -5.48, -5.10) CV[48,.] = (-5.88, -5.50, -5.20, -4.89) CV[49,.] = (-5.64, -5.33, -5.07, -4.77) CV[50,.] = (-5.60, -5.31, -5.05, -4.74) real colvector Tgrid Tgrid = (25\50\100\250\500) real scalar m1eff, m2eff, j, iT m1eff = max((min((m1,4)), 0)) m2eff = max((min((m2,2)), 1)) real scalar block block = (m2eff-1)*25 + m1eff*5 // pick T row by interpolation real rowvector cv cv = J(1,4,.) if (T <= 25) cv = CV[block+1, .] else if (T >= 500) cv = CV[block+5, .] else { real scalar lo, hi, w for (iT=1; iT<=4; iT++) { if (T >= Tgrid[iT] & T <= Tgrid[iT+1]) { lo = Tgrid[iT]; hi = Tgrid[iT+1] w = (T-lo)/(hi-lo) cv = (1-w)*CV[block+iT,.] + w*CV[block+iT+1,.] break } } } return(cv) } // -------------------------------------------------------------------- // Kongsted (2005) I(2)-to-I(1) LR transformation test // // Given an estimated I(2) VAR with rank (r, s1), and a user-supplied // p x q matrix G of candidate linear combinations, test H_0: sp(τ) = sp(G) // where τ = (β, β_⊥1). LR is asymptotically χ²(df) with // df = (p - r - s1) * q by Johansen (2006). // // This routine computes the LR by comparing the unrestricted I(2) log- // likelihood with the value obtained when τ is constrained to span G. // For practical purposes we implement the test as a Wald check: // statistic = T * trace( G⊥' β2 (β2' β2)^-1 β2' G⊥ ) // which is asymptotically χ²(df) under H_0 (Kurita 2011 §III). // -------------------------------------------------------------------- void _mixi12_kongsted(real matrix beta2, real matrix G, real scalar Tn, real scalar LR, real scalar df) { real matrix Gperp Gperp = _mixi12_orthcomp(G) real matrix M M = beta2' * Gperp LR = Tn * trace(M' * invsym(beta2'beta2) * M) df = (rows(Gperp) > 0 ? cols(Gperp) : 0) * cols(beta2) } // -------------------------------------------------------------------- // Stata-side wrappers (defined in the kernel so any sub-command can // call them after a single `_mixi12_mata` load). // -------------------------------------------------------------------- void _mixi12_haldrup_cv_runner(real scalar m1, real scalar m2, real scalar T) { real rowvector cv cv = _mixi12_haldrup_cv(m1, m2, T) st_numscalar("r(cv01)", cv[1]) st_numscalar("r(cv025)", cv[2]) st_numscalar("r(cv05)", cv[3]) st_numscalar("r(cv10)", cv[4]) } end