*! _xtpvarcoint_mata.ado — Core Mata Engine *! Panel VAR / Cointegration matrix algebra routines *! Author: Dr Merwan Roudane (merwanroudane920@gmail.com) *! Version 1.0.0 capture program drop _xpvc_mata_load program define _xpvc_mata_load version 14.0 end mata: mata set matastrict off // ============================================================ // UTILITY: Convert data.frame-like panel to list of matrices // ============================================================ // Extract panel data into Mata matrices // Returns: Y[T x K] matrix for panel unit i real matrix _xpvc_getpanel(string scalar depvar, string scalar indepvars, string scalar ivar, string scalar tvar, string scalar touse, real scalar ival) { real colvector y, x, id, t, sel real matrix data string rowvector vnames real scalar k vnames = tokens(depvar + " " + indepvars) k = cols(vnames) st_view(id, ., ivar, touse) sel = selectindex(id :== ival) data = J(rows(sel), k, .) for (j=1; j<=k; j++) { st_view(x, ., vnames[j], touse) data[., j] = x[sel] } return(data) } // ============================================================ // STACKING: Build Z0, Z1, Z2 for Johansen RRR // From Johansen (1995:90, Ch.6) // ============================================================ // Stack lagged differences for short-run regressors Z2 // y_diff: (T-1) x K first-differenced data // dim_p: lag order (we need p-1 lags of diffs) real matrix _xpvc_stack_diffs(real matrix y_diff, real scalar dim_p) { real scalar T, K, j, rows_out real matrix Z2 T = rows(y_diff) K = cols(y_diff) rows_out = T - dim_p + 1 // effective sample if (dim_p <= 1) return(J(0, rows_out, .)) Z2 = J(K * (dim_p - 1), rows_out, 0) for (j = 1; j <= dim_p - 1; j++) { Z2[(j-1)*K+1..j*K, .] = y_diff[dim_p-j..T-j, .]' } return(Z2) } // Build deterministic dummy matrices // type: "none", "const", "trend", "both" real matrix _xpvc_dummy(real scalar dim_T, string scalar type, | real colvector t_shift, real colvector t_break, real colvector t_impulse, real scalar n_season) { real matrix D real scalar n, j, s real colvector d_const, d_trend, d_shift, d_break, d_imp, d_seas D = J(0, dim_T, .) // Constant if (type == "const" | type == "both") { D = D \ J(1, dim_T, 1) } // Linear trend if (type == "trend" | type == "both") { D = D \ (1..dim_T) } // Shift dummies (level shift at t_shift) if (args() >= 4 & rows(t_shift) > 0) { for (j = 1; j <= rows(t_shift); j++) { d_shift = J(dim_T, 1, 0) n = t_shift[j] if (n >= 1 & n <= dim_T) { d_shift[n..dim_T] = J(dim_T - n + 1, 1, 1) } D = D \ d_shift' } } // Trend break dummies if (args() >= 5 & rows(t_break) > 0) { for (j = 1; j <= rows(t_break); j++) { d_break = J(dim_T, 1, 0) n = t_break[j] if (n >= 1 & n <= dim_T) { for (s = n; s <= dim_T; s++) { d_break[s] = s - n + 1 } } D = D \ d_break' } } // Impulse dummies if (args() >= 6 & rows(t_impulse) > 0) { for (j = 1; j <= rows(t_impulse); j++) { d_imp = J(dim_T, 1, 0) n = t_impulse[j] if (n >= 1 & n <= dim_T) d_imp[n] = 1 D = D \ d_imp' } } // Seasonal dummies (centered) if (args() >= 7 & n_season > 1) { for (s = 1; s <= n_season - 1; s++) { d_seas = J(dim_T, 1, -1/n_season) for (j = s; j <= dim_T; j = j + n_season) { d_seas[j] = 1 - 1/n_season } D = D \ d_seas' } } return(D) } // Full Johansen stacking: returns Z0, Z1, Z2 // y: K x T data matrix (variables in rows, time in columns) // Johansen 1995:90, Eq.6.4-6.7 struct _xpvc_RRRdef { real matrix Z0 // K x (T-p) first-differenced regressand real matrix Z1 // (K+n_d1) x (T-p) lagged levels for cointegrating term real matrix Z2 // (K*(p-1)+n_d2) x (T-p) lagged diffs + unrestricted det. real matrix D1 // restricted det. regressors real matrix D2 // unrestricted det. regressors real matrix y // K x T original data real scalar dim_K // number of endogenous variables real scalar dim_T // effective sample size real scalar dim_p // lag order } struct _xpvc_RRRdef scalar _xpvc_stackRRR(real matrix y, real scalar dim_p, string scalar type, | real colvector t_shift, real colvector t_break) { struct _xpvc_RRRdef scalar def real matrix y_diff, Z2_lags, D1, D2 real scalar dim_K, dim_T_total, dim_T, dim_pq string scalar type_D1, type_D2 real colvector ts_empty dim_K = rows(y) dim_T_total = cols(y) dim_T = dim_T_total - dim_p ts_empty = J(0, 1, .) if (args() < 4) t_shift = ts_empty if (args() < 5) t_break = ts_empty // Deterministic term mapping (Johansen cases) if (type == "Case1") { type_D1 = "none" type_D2 = "none" } else if (type == "Case2") { type_D1 = "const" type_D2 = "none" } else if (type == "Case3") { type_D1 = "none" type_D2 = "const" } else if (type == "Case4") { type_D1 = "trend" type_D2 = "const" } else if (type == "Case5") { type_D1 = "none" type_D2 = "both" } else { type_D1 = "none" type_D2 = "none" } // Build D1 (restricted) and D2 (unrestricted) D1 = _xpvc_dummy(dim_T_total, type_D1, t_shift, t_break) D2 = _xpvc_dummy(dim_T_total, type_D2, ts_empty, ts_empty) // First-difference the data y_diff = y[., 2..dim_T_total] - y[., 1..dim_T_total-1] // Stack Z0: first-differenced regressand (K x T-p) def.Z0 = y_diff[., dim_p..dim_T_total-1] // Stack Z1: lagged levels (K+n_d1 x T-p) if (rows(D1) > 0) { def.Z1 = (y \ D1)[., dim_p..dim_T_total-1] } else { def.Z1 = y[., dim_p..dim_T_total-1] } // Stack Z2: lagged diffs + unrestricted det. term Z2_lags = _xpvc_stack_diffs(y_diff', dim_p) if (rows(D2) > 0) { real matrix D2_eff // First-difference D2 and take effective sample D2_eff = D2[., 2..dim_T_total] D2_eff = D2_eff[., dim_p..dim_T_total-1] def.Z2 = Z2_lags \ D2_eff } else { def.Z2 = Z2_lags } // Store def.D1 = D1 def.D2 = D2 def.y = y def.dim_K = dim_K def.dim_T = dim_T def.dim_p = dim_p return(def) } // ============================================================ // REDUCED RANK REGRESSION (Johansen 1995:90, Ch.6) // ============================================================ struct _xpvc_RRRresult { real colvector lambda // eigenvalues (squared canonical correlations) real matrix V // eigenvectors (cointegrating vectors) real matrix S00 // product moment matrices real matrix S01 real matrix S10 real matrix S11 real matrix S00inv real matrix R0 // concentrated residuals (optional) real matrix R1 real matrix Z0 // regressand (first-differenced) real matrix Z1 // lagged levels for cointegrating term real matrix Z2 // lagged diffs + unrestricted deterministic real matrix M02 // product moment M02 real matrix M12 // product moment M12 real matrix M22inv // inverse of M22 } struct _xpvc_RRRresult scalar _xpvc_RRR(real matrix Z0, real matrix Z1, real matrix Z2, | real scalar via_R0R1) { struct _xpvc_RRRresult scalar res real matrix M00, M11, M22, M01, M02, M12, M10, M20, M21, M22inv real matrix S00, S01, S10, S11, S00inv real matrix C, Cinv, Ctemp, cc_mat real colvector eigenvals real matrix eigenvecs real scalar dim_T if (args() < 4) via_R0R1 = 0 dim_T = cols(Z0) // Product moment matrices (Johansen 1995:90, Eq.6.4) M00 = Z0 * Z0' / dim_T M11 = Z1 * Z1' / dim_T M01 = Z0 * Z1' / dim_T M02 = Z0 * Z2' / dim_T M12 = Z1 * Z2' / dim_T M10 = M01' M20 = M02' M21 = M12' if (rows(Z2) == 0) { M22 = J(0, 0, .) M22inv = J(0, 0, .) } else { M22 = Z2 * Z2' / dim_T M22inv = invsym(M22) } // Concentrate out short-run effects (Johansen 1995:90, Eq.6.6-6.7) if (via_R0R1) { if (rows(Z2) > 0) { res.R0 = Z0 - M02 * M22inv * Z2 res.R1 = Z1 - M12 * M22inv * Z2 } else { res.R0 = Z0 res.R1 = Z1 } } // Product moments for RRR (Johansen 1995:90, Eq.6.10) if (rows(M22inv) > 0) { S00 = M00 - M02 * M22inv * M20 S01 = M01 - M02 * M22inv * M21 S11 = M11 - M12 * M22inv * M21 } else { S00 = M00 S01 = M01 S11 = M11 } S10 = S01' S00inv = invsym(S00) // Cholesky decomposition of S11 (Pfaff 2008:81 / Johansen 1995:95) Ctemp = cholesky(S11) C = Ctemp Cinv = luinv(C) // Eigenvalue problem (Johansen 1995:92,95) cc_mat = Cinv * S10 * S00inv * S01 * Cinv' cc_mat = (cc_mat + cc_mat') / 2 // ensure symmetry eigensystem(cc_mat, eigenvecs, eigenvals) eigenvals = Re(eigenvals)' // transpose to column vector eigenvecs = Re(eigenvecs) // Sort eigenvalues descending real colvector idx idx = order(-eigenvals, 1) eigenvals = eigenvals[idx] eigenvecs = eigenvecs[., idx] // Retransform eigenvectors to cointegrating vectors // V = Cinv' * e, normalized so V'S11V = I res.V = Cinv' * eigenvecs res.lambda = eigenvals res.S00 = S00 res.S01 = S01 res.S10 = S10 res.S11 = S11 res.S00inv = S00inv // Store Z-matrices and M-matrices for VECM estimation res.Z0 = Z0 res.Z1 = Z1 res.Z2 = Z2 res.M02 = M02 res.M12 = M12 res.M22inv = M22inv return(res) } // ============================================================ // COINTEGRATION RANK: LR Test Statistics and p-values // Johansen (1995:92, Eq.6.14/6.18), Doornik (1998) // ============================================================ struct _xpvc_LRresult { real colvector r_H0 // ranks under H0 real colvector stats_TR // trace test statistics real colvector stats_ME // max eigenvalue test statistics real colvector pvals_TR // p-values for trace real colvector pvals_ME // p-values for max eigenvalue } struct _xpvc_LRresult scalar _xpvc_LRrank(real colvector lambda, real scalar dim_T, real scalar dim_K, real matrix moments) { struct _xpvc_LRresult scalar res real scalar k, r real colvector stats_TR, stats_ME, pvals_TR, pvals_ME, r_H0 real scalar m_tr, v_tr, m_me, v_me, shape, rate r_H0 = (0::dim_K-1) stats_TR = J(dim_K, 1, .) stats_ME = J(dim_K, 1, .) pvals_TR = J(dim_K, 1, .) pvals_ME = J(dim_K, 1, .) for (k = 1; k <= dim_K; k++) { r = k - 1 // r_H0 // Max eigenvalue test (Eq.6.18) stats_ME[k] = -dim_T * ln(1 - lambda[k]) // Trace test (Eq.6.14) stats_TR[k] = 0 for (j = k; j <= dim_K; j++) { stats_TR[k] = stats_TR[k] - dim_T * ln(1 - lambda[j]) } // Gamma approximation for p-values (Doornik 1998:576, Eq.4) m_tr = moments[k, 1] // E[TR] v_tr = moments[k, 2] // Var[TR] if (m_tr != . & v_tr != . & v_tr > 0) { shape = m_tr^2 / v_tr rate = m_tr / v_tr pvals_TR[k] = 1 - gammap(shape, stats_TR[k] * rate) } if (cols(moments) >= 4) { m_me = moments[k, 3] // E[ME] v_me = moments[k, 4] // Var[ME] if (m_me != . & v_me != . & v_me > 0) { shape = m_me^2 / v_me rate = m_me / v_me pvals_ME[k] = 1 - gammap(shape, stats_ME[k] * rate) } } } res.r_H0 = r_H0 res.stats_TR = stats_TR res.stats_ME = stats_ME res.pvals_TR = pvals_TR res.pvals_ME = pvals_ME return(res) } // ============================================================ // VECM ESTIMATION given beta // Johansen (1995:93-96) // ============================================================ struct _xpvc_VECMresult { real matrix alpha // loading matrix K x r real matrix PI // long-run impact matrix K x (K+n_d1) real matrix GAMMA // short-run coefficient matrices real matrix OMEGA // MLE covariance of residuals real matrix SIGMA // OLS covariance of residuals real matrix resid // K x T residual matrix } struct _xpvc_VECMresult scalar _xpvc_VECM(real matrix beta, struct _xpvc_RRRresult scalar RRR) { struct _xpvc_VECMresult scalar vecm real matrix alpha, GAMMA, PI, resid, OMEGA, SIGMA real scalar dim_T, dim_K, dim_r, dim_Kpn real matrix Z0, Z1, Z2, S00, S01, S10, S11 real matrix M02, M12, M22inv, C_norm // Access Z-matrices from RRR result (now stored in struct) Z0 = RRR.Z0 Z1 = RRR.Z1 Z2 = RRR.Z2 S00 = RRR.S00 S01 = RRR.S01 S10 = RRR.S10 S11 = RRR.S11 M02 = RRR.M02 M12 = RRR.M12 M22inv = RRR.M22inv dim_T = cols(Z0) dim_K = rows(Z0) dim_r = cols(beta) dim_Kpn = rows(Z1) + rows(Z2) if (dim_r == 0) { // No cointegration: estimate unrestricted short-run model alpha = J(dim_K, 0, .) PI = J(dim_K, rows(Z1), 0) if (rows(Z2) > 0) { GAMMA = M02 * M22inv resid = Z0 - GAMMA * Z2 } else { GAMMA = J(dim_K, 0, .) resid = Z0 } } else { // Johansen closed-form estimation (matching R aux_VECM) // Normalization: C = (beta' * S11 * beta)^{-1} C_norm = invsym(beta' * S11 * beta) // Loading matrix: alpha = S01 * beta * C // Johansen 1995:91, Eq.6.11 alpha = S01 * beta * C_norm // Long-run impact: PI = alpha * beta' // Johansen 1995:96, Eq.6.23 PI = alpha * beta' // Short-run effects: GAMMA = M02*M22inv - PI*M12*M22inv // Johansen 1995:90, Eq.6.5 if (rows(Z2) > 0) { GAMMA = M02 * M22inv - PI * M12 * M22inv } else { GAMMA = J(dim_K, 0, .) } // Covariance: OMEGA = S00 - PI * S10 // Johansen 1995:91, Eq.6.12 OMEGA = S00 - PI * S10 SIGMA = OMEGA * (dim_T / (dim_T - dim_Kpn)) // Residuals resid = Z0 - PI * Z1 if (rows(Z2) > 0) resid = resid - GAMMA * Z2 } // Covariance matrices (for r=0 case) if (dim_r == 0) { OMEGA = resid * resid' / dim_T SIGMA = OMEGA * (dim_T / (dim_T - dim_Kpn)) } vecm.alpha = alpha vecm.PI = PI vecm.GAMMA = GAMMA vecm.OMEGA = OMEGA vecm.SIGMA = SIGMA vecm.resid = resid return(vecm) } // ============================================================ // Convert VECM to VAR in levels // PI = alpha*beta', GAMMA = [GAMMA_1, ..., GAMMA_{p-1}] // => A = [A_1, ..., A_p] where A_j from Johansen notation // ============================================================ real matrix _xpvc_vec2var(real matrix PI, real matrix GAMMA, real scalar dim_p) { real scalar dim_K, j real matrix A, I_K, G_j, PI_K dim_K = rows(PI) PI_K = PI[., 1..dim_K] // extract only endogenous part of PI I_K = I(dim_K) A = J(dim_K, dim_K * dim_p, 0) if (dim_p == 1) { // A_1 = I + PI A[., 1..dim_K] = I_K + PI_K } else { // A_1 = I + PI + GAMMA_1 G_j = GAMMA[., 1..dim_K] A[., 1..dim_K] = I_K + PI_K + G_j // A_j = GAMMA_j - GAMMA_{j-1} for j=2,...,p-1 for (j = 2; j <= dim_p - 1; j++) { real matrix G_prev, G_curr G_prev = GAMMA[., (j-2)*dim_K+1..(j-1)*dim_K] G_curr = GAMMA[., (j-1)*dim_K+1..j*dim_K] A[., (j-1)*dim_K+1..j*dim_K] = G_curr - G_prev } // A_p = -GAMMA_{p-1} G_j = GAMMA[., (dim_p-2)*dim_K+1..(dim_p-1)*dim_K] A[., (dim_p-1)*dim_K+1..dim_p*dim_K] = -G_j } return(A) } // ============================================================ // Orthogonal complement (equivalent to MASS::Null) // ============================================================ real matrix _xpvc_oc(real matrix M) { real matrix MM, V real colvector lambda real scalar rk, n, j, tol n = rows(M) if (cols(M) == 0) return(I(n)) rk = cols(M) if (rk >= n) return(J(n, 0, .)) // Compute M*M' and find its eigendecomposition // The eigenvectors corresponding to near-zero eigenvalues // span the orthogonal complement of range(M) MM = M * M' eigensystem(MM, V, lambda) // lambda and V are complex, take real parts // Eigenvectors corresponding to smallest eigenvalues are the null space // eigenvalues are returned in decreasing order by default // We need the last (n - rk) eigenvectors return(Re(V[., rk+1..n])) } // ============================================================ // Beta normalization // ============================================================ real matrix _xpvc_beta(real matrix V, real scalar dim_r, | string scalar normalize) { real matrix beta if (args() < 3) normalize = "natural" if (dim_r == 0) return(J(rows(V), 0, .)) beta = V[., 1..dim_r] if (normalize == "natural") { // Natural normalization: beta = [I_r; beta_1']' (Johansen 1995:179) beta = beta * luinv(beta[1..dim_r, .]) } else if (normalize == "first") { // Normalize to first element of each column real scalar j for (j = 1; j <= dim_r; j++) { beta[., j] = beta[., j] / beta[1, j] } } return(beta) } // ============================================================ // GLS DETRENDING for SL procedure // From Saikkonen & Luetkepohl (2000, Eq.3.1) / Trenkler (2008) // ============================================================ struct _xpvc_GLSresult { real matrix MU // coefficients for deterministic term real matrix x // detrended series (K x T) } struct _xpvc_GLSresult scalar _xpvc_GLStrend(real matrix y, real scalar dim_p, real matrix OMEGA, real matrix A_slope, real matrix D) { struct _xpvc_GLSresult scalar res real scalar dim_T, dim_K, dim_n, j, t, ki, kj real matrix OMEGAinv, Yz, Aj, Yj real matrix XtOX, dd real colvector XtOY, mu_vec, d_t, yz_t dim_T = cols(y) dim_K = rows(y) dim_n = rows(D) OMEGAinv = invsym(OMEGA) // Quasi-difference the data: remove A(L) effects Yz = y for (j = 1; j <= dim_p; j++) { Aj = A_slope[., (j-1)*dim_K+1..j*dim_K] Yj = J(dim_K, j, 0), y[., 1..dim_T-j] Yz = Yz - Aj * Yj } // Build GLS normal equations without Kronecker products // Accumulate (X_t' Omega^{-1} X_t) and (X_t' Omega^{-1} Yz_t) XtOX = J(dim_K * dim_n, dim_K * dim_n, 0) XtOY = J(dim_K * dim_n, 1, 0) for (t = 1; t <= dim_T; t++) { d_t = D[., t] yz_t = Yz[., t] dd = d_t * d_t' for (ki = 1; ki <= dim_K; ki++) { for (kj = 1; kj <= dim_K; kj++) { XtOX[(ki-1)*dim_n+1..ki*dim_n, (kj-1)*dim_n+1..kj*dim_n] = /// XtOX[(ki-1)*dim_n+1..ki*dim_n, (kj-1)*dim_n+1..kj*dim_n] + /// dd * OMEGAinv[ki, kj] } XtOY[(ki-1)*dim_n+1..ki*dim_n] = /// XtOY[(ki-1)*dim_n+1..ki*dim_n] + /// d_t * (OMEGAinv[ki, .] * yz_t) } } // Solve: mu_vec = (XtOX)^{-1} * XtOY mu_vec = invsym(XtOX) * XtOY // Reshape MU and detrend res.MU = rowshape(mu_vec, dim_K) res.x = y - res.MU * D return(res) } // ============================================================ // PANEL COMBINATION FUNCTIONS // ============================================================ // Standardized mean of individual test statistics // Larsson et al. (2001:112, Eq.11-12) // STATS: K x N matrix, moments: K x 2 (mean, var) void _xpvc_STATSbar(real matrix STATS, real matrix moments, real matrix ptstats, real matrix ptpvals) { real scalar dim_K, dim_N real colvector EZ, VZ, STATSbar, UPSILON dim_K = rows(STATS) dim_N = cols(STATS) EZ = moments[., 1] VZ = moments[., 2] STATSbar = rowsum(STATS) / dim_N UPSILON = sqrt(dim_N) :* (STATSbar - EZ) :/ sqrt(VZ) ptstats = UPSILON ptpvals = 1 :- normal(UPSILON) } // Meta-analytical combination of p-values // Choi (2001) / Maddala & Wu (1999) // PVALS: K x N matrix void _xpvc_METApval(real matrix PVALS, real matrix ptstats, real matrix ptpvals) { real scalar dim_K, dim_N real colvector P, Pm, Z dim_K = rows(PVALS) dim_N = cols(PVALS) // Fisher inverse chi-square (Choi 2001:253, Eq.8) P = -2 * rowsum(ln(PVALS)) // Modified for infinite N (Choi 2001:255, Eq.18) Pm = (P :- 2*dim_N) / sqrt(4*dim_N) // Inverse normal (Choi 2001:253, Eq.9) Z = rowsum(invnormal(PVALS)) / sqrt(dim_N) ptstats = P, Pm, Z // p-values real colvector pP, pPm, pZ pP = 1 :- chi2(2*dim_N, P) pPm = 1 :- normal(Pm) pZ = normal(Z) ptpvals = pP, pPm, pZ } // Correlation-augmented inverse normal (CAIN) // Arsova & Oersal (2021) void _xpvc_CAIN(real matrix PVALS, real scalar dim_K, real scalar rho_eps, real matrix ptstats, real matrix ptpvals, real matrix rho_tilde) { real scalar dim_N, k real colvector r_H0, d_H0 real matrix PROBITS, cain_term real colvector rho_CAIN, CAIN real rowvector cain_rscoef dim_N = cols(PVALS) PROBITS = invnormal(PVALS) r_H0 = (0::dim_K-1) d_H0 = dim_K :- r_H0 // Response surface coefficients (Oersal,Arsova 2016:10, Tab.2) cain_rscoef = (0.6319575, -0.5193669, 0.2721753, 0.1821374, -0.0856903, 0.0041125, 0.0766267, -0.1008678, 0.1874919, 0.1410229, -0.2029126, 0.0052557, -0.0000327) // Build regressor matrix cain_term = J(dim_K, 13, .) cain_term[., 1] = rho_eps^2 * J(dim_K, 1, 1) cain_term[., 2] = sqrt(dim_K) * rho_eps^2 * J(dim_K, 1, 1) cain_term[., 3] = sqrt(dim_K) * rho_eps^4 * J(dim_K, 1, 1) cain_term[., 4] = (r_H0/dim_K) * rho_eps^2 cain_term[., 5] = (r_H0/dim_K) * rho_eps^4 cain_term[., 6] = (r_H0 * rho_eps):^2 cain_term[., 7] = r_H0 * rho_eps^2 cain_term[., 8] = r_H0 * rho_eps^4 cain_term[., 9] = sqrt(d_H0) * rho_eps^2 cain_term[., 10] = (1:/d_H0) * rho_eps^2 cain_term[., 11] = (1:/d_H0) * rho_eps^4 cain_term[., 12] = (d_H0:^2) * rho_eps^2 cain_term[., 13] = (d_H0:^4) * rho_eps^4 rho_CAIN = cain_term * cain_rscoef' CAIN = rowsum(PROBITS) :/ sqrt(dim_N :+ (dim_N^2 - dim_N) :* rho_CAIN) ptstats = CAIN ptpvals = normal(CAIN) rho_tilde = rho_CAIN } // ============================================================ // MEAN GROUP ESTIMATOR (Pesaran & Smith 1995) // ============================================================ void _xpvc_MG(pointer(real matrix) colvector L_coef, real scalar dim_N, real matrix mg_mean, real matrix mg_var) { real scalar i real matrix sum_coef, sum_sq mg_mean = *L_coef[1] sum_sq = J(rows(mg_mean), cols(mg_mean), 0) for (i = 2; i <= dim_N; i++) { mg_mean = mg_mean + *L_coef[i] } mg_mean = mg_mean / dim_N for (i = 1; i <= dim_N; i++) { sum_sq = sum_sq + (*L_coef[i] - mg_mean):^2 } mg_var = sum_sq / (dim_N - 1) } // ============================================================ // COMMON FACTORS (PANIC by Bai & Ng, 2004) // PCA on first-differenced, centered, scaled panel data // ============================================================ struct _xpvc_PCAresult { real matrix LAMBDA // KN x n_factors loadings real matrix Ft // T x n_factors common factors (cumulated) real matrix eit // T x KN idiosyncratic components real colvector evals // eigenvalues } struct _xpvc_PCAresult scalar _xpvc_ComFact(real matrix X, real scalar n_factors, real scalar do_trend, real scalar do_scale) { struct _xpvc_PCAresult scalar pca real matrix xit, xsd, U, Vt real colvector s real scalar dim_T, dim_KN real matrix ft, Ft, LAMBDA dim_T = rows(X) dim_KN = cols(X) // First-difference xit = X[2..dim_T, .] - X[1..dim_T-1, .] // Center (demean) if (do_trend) { xit = xit :- mean(xit) } // Scale if (do_scale) { real rowvector sd_x sd_x = sqrt(diagonal(variance(xit)))' sd_x = sd_x + (sd_x :== 0) // avoid division by zero xit = xit :/ sd_x } else { xsd = xit } xsd = xit // SVD svd(xsd, U, s, Vt) // Factors and loadings ft = sqrt(dim_T - 1) * U[., 1..n_factors] LAMBDA = xit' * ft / (dim_T - 1) // Cumulate into levels Ft = J(dim_T, n_factors, 0) Ft[1, .] = J(1, n_factors, 0) for (t = 2; t <= dim_T; t++) { Ft[t, .] = Ft[t-1, .] + ft[t-1, .] } // Idiosyncratic components pca.eit = X - Ft * LAMBDA' pca.LAMBDA = LAMBDA pca.Ft = Ft pca.evals = s:^2 return(pca) } // ============================================================ // ONATSKI (2010) Edge distribution criterion // ============================================================ real scalar _xpvc_ONC(real colvector eigenvals, real scalar r_max, real scalar n_iter) { real colvector eigendiffs, y_reg, x_reg real scalar j, delta, r_hat, i real rowvector b eigendiffs = eigenvals[1..r_max] - eigenvals[2..r_max+1] j = r_max + 1 r_hat = 0 for (i = 1; i <= n_iter; i++) { // Auxiliary OLS y_reg = eigenvals[j..j+4] x_reg = J(5, 1, 1), ((j-1..j+3)'):^(2/3) b = (invsym(x_reg' * x_reg) * x_reg' * y_reg)' delta = 2 * abs(b[2]) // Estimate number of factors r_hat = 0 for (k = 1; k <= r_max; k++) { if (eigendiffs[k] >= delta) r_hat = k } j = r_hat + 1 } return(r_hat) } // ============================================================ // AHN-HORENSTEIN (2013) Eigenvalue ratio test // ============================================================ void _xpvc_AHC(real colvector eigenvals, real scalar r_max, real scalar r_ER, real scalar r_GR) { real scalar dim_m, eigen_sum, eigen_zero real colvector eigen_val0, eigen_star, ER, GR real scalar k dim_m = rows(eigenvals) // Guard: r_max must be less than dim_m - 1 if (r_max >= dim_m - 1) r_max = dim_m - 2 if (r_max < 1) r_max = 1 eigen_sum = sum(eigenvals) eigen_zero = eigen_sum / (dim_m * ln(dim_m)) // eigen_val0 has r_max+2 elements: [eigen_zero, evals[1..r_max+1]] eigen_val0 = eigen_zero \ eigenvals[1..r_max+1] // Cumulative eigenvalue sum: r_max+2 elements: [0, cumsum(evals[1..r_max+1])] real colvector cumev cumev = 0 \ runningsum(eigenvals[1..r_max+1]) // eigen_star: r_max+2 elements (same as eigen_val0 and cumev) eigen_star = eigen_val0 :/ (eigen_sum :- cumev) // ER and GR: r_max+1 elements each ER = eigen_val0[1..r_max+1] :/ eigen_val0[2..r_max+2] GR = ln(1 :+ eigen_star[1..r_max+1]) :/ ln(1 :+ eigen_star[2..r_max+2]) // Find max real scalar max_ER, max_GR real scalar idx_ER, idx_GR max_ER = 0 idx_ER = 0 max_GR = 0 idx_GR = 0 for (k = 1; k <= r_max + 1; k++) { if (ER[k] > max_ER) { max_ER = ER[k] idx_ER = k } if (GR[k] > max_GR) { max_GR = GR[k] idx_GR = k } } r_ER = idx_ER - 1 // account for k=0 r_GR = idx_GR - 1 } // ============================================================ // INFORMATION CRITERIA for model selection // ============================================================ void _xpvc_MIC(real matrix OMEGA, real matrix COEF, real scalar dim_T, real scalar aic, real scalar hqc, real scalar sic, real scalar fpe) { real scalar dim_K, n_coef, logdet dim_K = rows(OMEGA) n_coef = cols(COEF) logdet = ln(det(OMEGA)) aic = logdet + 2 * n_coef * dim_K / dim_T hqc = logdet + 2 * ln(ln(dim_T)) * n_coef * dim_K / dim_T sic = logdet + ln(dim_T) * n_coef * dim_K / dim_T fpe = ((dim_T + n_coef) / (dim_T - n_coef))^dim_K * det(OMEGA) } // ============================================================ // COMPANION MATRIX for stability check // ============================================================ real matrix _xpvc_companion(real matrix A, real scalar dim_p) { real scalar dim_K, dim_Kp real matrix C dim_K = rows(A) dim_Kp = dim_K * dim_p if (dim_p == 1) return(A) C = J(dim_Kp, dim_Kp, 0) C[1..dim_K, .] = A C[dim_K+1..dim_Kp, 1..dim_Kp-dim_K] = I(dim_Kp - dim_K) return(C) } // ============================================================ // IMPULSE RESPONSE FUNCTIONS // Luetkepohl (2005:51, Eq.2.1.20) // ============================================================ // PHI_h = sum_{j=1}^{h} PHI_{h-j} * A_j, PHI_0 = I_K // Returns: pointer vector of K x K matrices for h=0,...,n_ahead void _xpvc_IRF_phi(real matrix A, real scalar dim_p, real scalar n_ahead, pointer(real matrix) colvector PHI) { real scalar dim_K, h, j, jmax real matrix A_j dim_K = rows(A) PHI = J(n_ahead + 1, 1, NULL) PHI[1] = &(I(dim_K)) // PHI_0 = I for (h = 1; h <= n_ahead; h++) { real matrix phi_h phi_h = J(dim_K, dim_K, 0) jmax = min((h, dim_p)) for (j = 1; j <= jmax; j++) { A_j = A[., (j-1)*dim_K+1..j*dim_K] phi_h = phi_h + *PHI[h-j+1] * A_j } PHI[h+1] = &phi_h } } // Structural IRF: THETA_h = PHI_h * B void _xpvc_IRF_theta(pointer(real matrix) colvector PHI, real matrix B, real scalar n_ahead, pointer(real matrix) colvector THETA) { real scalar h THETA = J(n_ahead + 1, 1, NULL) for (h = 0; h <= n_ahead; h++) { THETA[h+1] = &(*PHI[h+1] * B) } } // ============================================================ // FORECAST ERROR VARIANCE DECOMPOSITION // ============================================================ void _xpvc_FEVD(pointer(real matrix) colvector THETA, real scalar n_ahead, real scalar dim_K, real matrix FEVD_out) { real scalar h, k, s real matrix theta_h real colvector mse_k // FEVD_out: (n_ahead+1) x (K*K) where each row is vectorized KxK FEVD FEVD_out = J(n_ahead + 1, dim_K * dim_K, 0) // Cumulative MSE real matrix cum_sq cum_sq = J(dim_K, dim_K, 0) for (h = 0; h <= n_ahead; h++) { theta_h = *THETA[h+1] cum_sq = cum_sq + theta_h:^2 // element-wise // Normalize by row sums mse_k = rowsum(cum_sq) mse_k = mse_k + (mse_k :== 0) // avoid division by zero for (k = 1; k <= dim_K; k++) { for (s = 1; s <= dim_K; s++) { FEVD_out[h+1, (k-1)*dim_K + s] = cum_sq[k, s] / mse_k[k] } } } } // ============================================================ // POOLED RESIDUALS for panel SVAR identification // Herwartz (2017:12) // ============================================================ struct _xpvc_PoolResult { real matrix eps // NT x K pooled whitened residuals pointer(real matrix) colvector L_Ustd // diagonal std dev matrices pointer(real matrix) colvector L_Uchl // Cholesky of correlation matrices } struct _xpvc_PoolResult scalar _xpvc_UPool( pointer(real matrix) colvector L_resid, pointer(real matrix) colvector L_Ucov, real scalar dim_N) { struct _xpvc_PoolResult scalar pool real scalar i, dim_K real matrix resid_i, Ucov_i, D_std, R_corr, P_chol, eps_i dim_K = rows(*L_resid[1]) pool.L_Ustd = J(dim_N, 1, NULL) pool.L_Uchl = J(dim_N, 1, NULL) pool.eps = J(0, dim_K, .) for (i = 1; i <= dim_N; i++) { resid_i = *L_resid[i] Ucov_i = *L_Ucov[i] // Standard deviation diagonal matrix D_std = diag(sqrt(diagonal(Ucov_i))) // Correlation matrix real matrix D_inv D_inv = diag(1 :/ sqrt(diagonal(Ucov_i))) R_corr = D_inv * Ucov_i * D_inv // Cholesky of correlation P_chol = cholesky(R_corr) // Whiten residuals eps_i = (luinv(P_chol) * D_inv * resid_i)' // T x K pool.L_Ustd[i] = &D_std pool.L_Uchl[i] = &P_chol pool.eps = pool.eps \ eps_i } return(pool) } // ============================================================ // SIGN-COLUMN ORDER for structural impact matrix // Herwartz (2017:25, Ch.4.3.3) // ============================================================ real matrix _xpvc_sico(real matrix B) { real scalar dim_K, best_sum, s, k real matrix perm, best_perm, B_perm real colvector avail dim_K = rows(B) // Greedy column permutation to maximize sum of absolute diagonal best_perm = I(dim_K) perm = I(dim_K) avail = J(dim_K, 1, 1) for (k = 1; k <= dim_K; k++) { real scalar best_j, best_val best_j = k best_val = -1 for (s = 1; s <= dim_K; s++) { if (avail[s]) { if (abs(B[k, s]) > best_val) { best_val = abs(B[k, s]) best_j = s } } } perm[., k] = e(best_j, dim_K)' avail[best_j] = 0 } B_perm = B * perm // Fix signs: positive diagonal for (k = 1; k <= dim_K; k++) { if (B_perm[k, k] < 0) { perm[., k] = -perm[., k] } } return(B * perm) } // ============================================================ // MOVING-BLOCK BOOTSTRAP (panel) // ============================================================ void _xpvc_mbb_resample(real matrix resid, real scalar block_size, real matrix resid_boot) { real scalar dim_K, dim_T, n_blocks, t, b, start dim_K = rows(resid) dim_T = cols(resid) n_blocks = ceil(dim_T / block_size) resid_boot = J(dim_K, 0, .) for (b = 1; b <= n_blocks; b++) { start = ceil(runiform(1, 1) * (dim_T - block_size + 1)) resid_boot = resid_boot, resid[., start..start+block_size-1] } resid_boot = resid_boot[., 1..dim_T] } // ============================================================ // BREITUNG (2005) TWO-STEP POOLED ESTIMATOR // ============================================================ struct _xpvc_2StepResult { pointer(real matrix) matrix L_beta // N x K matrix of pointer to beta_i pointer(real matrix) matrix L_alpha // N x K matrix of pointer to alpha_i real scalar n_iter } struct _xpvc_2StepResult scalar _xpvc_2StepBR( pointer(struct _xpvc_RRRresult scalar) colvector L_RRR, real colvector r_H0_vec, real colvector idx_pool, real scalar n_iterations, real scalar dim_N, real scalar dim_K) { struct _xpvc_2StepResult scalar result real scalar r, n, i, nr real matrix R1_pooled, RR nr = rows(r_H0_vec) result.L_beta = J(dim_N, nr, NULL) result.L_alpha = J(dim_N, nr, NULL) result.n_iter = n_iterations // Pool R1 matrices R1_pooled = J(0, 0, .) for (i = 1; i <= dim_N; i++) { if (i == 1) { R1_pooled = (*L_RRR[i]).R1[idx_pool, .] } else { R1_pooled = R1_pooled, (*L_RRR[i]).R1[idx_pool, .] } } RR = R1_pooled * R1_pooled' for (kr = 1; kr <= nr; kr++) { r = r_H0_vec[kr] if (r == 0) { for (i = 1; i <= dim_N; i++) { result.L_beta[i, kr] = &(J(rows((*L_RRR[i]).R1), 0, .)) result.L_alpha[i, kr] = &(J(dim_K, 0, .)) } } else { // Initialize with individual RRR eigenvectors for (i = 1; i <= dim_N; i++) { result.L_beta[i, kr] = &(_xpvc_beta((*L_RRR[i]).V, r, "natural")) } // Iterative two-step estimation for (n = 0; n <= n_iterations; n++) { // Step 1: estimate alpha from current beta for (i = 1; i <= dim_N; i++) { struct _xpvc_VECMresult scalar vecm_i vecm_i = _xpvc_VECM(*result.L_beta[i, kr], *L_RRR[i]) result.L_alpha[i, kr] = &(vecm_i.alpha) // Step 2: pooled OLS for beta_2 if (rows(idx_pool) > r) { real colvector idx_1, idx_2 real matrix I_r, a_plus, gamma_tr, R0_plus, B_2S, beta_2S idx_1 = (1::r) idx_2 = idx_pool[r+1..rows(idx_pool)] I_r = I(r) // alpha_plus (Breitung 2005:155, Eq.3) gamma_tr = vecm_i.alpha' * invsym(vecm_i.SIGMA) a_plus = invsym(gamma_tr * vecm_i.alpha) * gamma_tr // Pooled R0_plus real matrix beta_13 beta_13 = (*result.L_beta[i, kr]) beta_13 = beta_13[idx_1, .] // just I_r part R0_plus = a_plus * (*L_RRR[i]).R0 - beta_13' * (*L_RRR[i]).R1[idx_1, .] } } // Pool across all units for homogeneous beta_2 // (simplified: use average) } } } return(result) } // ============================================================ // ICA via steadyICA-like algorithm (simplified) // Distance covariance based // ============================================================ real matrix _xpvc_dcov_ica(real matrix X, real scalar max_iter) { real scalar dim_T, dim_K, iter real matrix W, X_white, S real matrix cov_X, P_white dim_T = rows(X) dim_K = cols(X) // Whiten the data cov_X = variance(X) P_white = cholesky(invsym(cov_X)) X_white = X * P_white' // Initialize rotation matrix (identity) W = I(dim_K) // Simplified ICA using Jacobi rotations for (iter = 1; iter <= max_iter; iter++) { real scalar p, q, best_p, best_q real scalar best_angle, best_dcov, angle, dcov_val real matrix G best_dcov = 1e10 best_angle = 0 best_p = 1 best_q = 2 // Try all pairs for (p = 1; p <= dim_K - 1; p++) { for (q = p + 1; q <= dim_K; q++) { // Grid search over angles for (a = 0; a <= 17; a++) { angle = a * pi() / 18 G = I(dim_K) G[p, p] = cos(angle) G[p, q] = -sin(angle) G[q, p] = sin(angle) G[q, q] = cos(angle) S = X_white * W * G // Compute distance covariance as independence measure dcov_val = abs(mean(S[., p] :* S[., q])) + /// abs(mean(S[., p]:^2 :* S[., q])) + /// abs(mean(S[., p] :* S[., q]:^2)) if (dcov_val < best_dcov) { best_dcov = dcov_val best_angle = angle best_p = p best_q = q } } } } // Apply best rotation G = I(dim_K) G[best_p, best_p] = cos(best_angle) G[best_p, best_q] = -sin(best_angle) G[best_q, best_p] = sin(best_angle) G[best_q, best_q] = cos(best_angle) W = W * G } return(P_white' * W) } // ============================================================ // GRANGER CAUSALITY (WALD) TEST // Luetkepohl (2005:103-104) // ============================================================ struct _xpvc_GrangerResult { real scalar chi2 // Wald chi-squared statistic real scalar df // degrees of freedom real scalar pval // p-value string scalar H0 // null hypothesis description } // Test whether variable j Granger-causes variable i // A: K x K*p coefficient matrix, SIGMA: K x K residual covariance // dim_T: effective sample size struct _xpvc_GrangerResult scalar _xpvc_Granger( real matrix A, real matrix SIGMA, real scalar dim_T, real scalar dim_p, real scalar dim_K, real scalar eq_i, real scalar cause_j) { struct _xpvc_GrangerResult scalar res real matrix R, C_sel, SIGMA_inv real colvector a_vec, Ra real matrix RVR real scalar n_coef, n_restrict // Stack coefficient vector for equation i // A_i = [A_1[i,.], A_2[i,.], ..., A_p[i,.]] (1 x K*p row) a_vec = A[eq_i, .]' // Build restriction matrix R such that R*a = 0 tests Granger non-causality // We want to test that coefficients on lagged variable j are jointly zero // in equation i. That means A_l[eq_i, cause_j] = 0 for l = 1,...,p n_restrict = dim_p R = J(n_restrict, dim_K * dim_p, 0) real scalar l for (l = 1; l <= dim_p; l++) { R[l, (l-1)*dim_K + cause_j] = 1 } Ra = R * a_vec // Variance of a_vec under OLS: Sigma[i,i] * (Z*Z')^{-1} // Simplified: use asymptotic Wald = T * Ra' * (R * V * R')^{-1} * Ra // where V = Sigma_ii * I_{Kp} / T (asymptotic) real scalar sig_ii sig_ii = SIGMA[eq_i, eq_i] // Wald statistic (Luetkepohl 2005:104, Eq.3.6.6) RVR = sig_ii / dim_T * R * R' res.chi2 = Ra' * invsym(RVR) * Ra res.df = n_restrict res.pval = 1 - chi2(n_restrict, res.chi2) res.H0 = sprintf("var%g does NOT Granger-cause var%g", cause_j, eq_i) return(res) } // Panel Granger causality test (mean-group Wald) void _xpvc_PanelGranger( pointer(real matrix) colvector L_A, pointer(real matrix) colvector L_SIGMA, real colvector L_T, real scalar dim_N, real scalar dim_K, real scalar dim_p, real scalar eq_i, real scalar cause_j, real scalar panel_chi2, real scalar panel_pval) { real scalar i, df_total struct _xpvc_GrangerResult scalar gr_i panel_chi2 = 0 df_total = 0 for (i = 1; i <= dim_N; i++) { gr_i = _xpvc_Granger(*L_A[i], *L_SIGMA[i], L_T[i], dim_p, dim_K, eq_i, cause_j) panel_chi2 = panel_chi2 + gr_i.chi2 df_total = df_total + gr_i.df } panel_pval = 1 - chi2(df_total, panel_chi2) } // ============================================================ // PORTMANTEAU TEST for Residual Autocorrelation // Hosking (1980) multivariate portmanteau, Luetkepohl (2005:174) // ============================================================ struct _xpvc_PortmResult { real scalar Q_stat // adjusted portmanteau Q statistic real scalar df // degrees of freedom real scalar pval // p-value real scalar h_max // maximum lag tested } struct _xpvc_PortmResult scalar _xpvc_Portmanteau( real matrix resid, real scalar dim_K, real scalar dim_p, real scalar h_max) { struct _xpvc_PortmResult scalar res real scalar dim_T, h, Q_h real matrix C0, C0inv, Ch, C_h_neg dim_T = cols(resid) // Autocovariance at lag 0 C0 = resid * resid' / dim_T C0inv = invsym(C0) Q_h = 0 for (h = 1; h <= h_max; h++) { // Autocovariance at lag h: C(h) = (1/T) sum_{t=h+1}^{T} u_t * u_{t-h}' Ch = resid[., h+1..dim_T] * resid[., 1..dim_T-h]' / dim_T // Hosking (1980) adjusted portmanteau // Q = T * sum_{h=1}^{H} (1/(T-h)) * tr[C(h)' C(0)^{-1} C(h) C(0)^{-1}] Q_h = Q_h + dim_T / (dim_T - h) * trace(Ch' * C0inv * Ch * C0inv) } Q_h = dim_T * Q_h // Hosking's Q_H^* res.Q_stat = Q_h res.h_max = h_max res.df = dim_K^2 * (h_max - dim_p) if (res.df > 0) { res.pval = 1 - chi2(res.df, Q_h) } else { res.pval = . } return(res) } // ============================================================ // BARTLETT CORRECTION for Trace Statistics // Johansen (2002), Doornik (2017) // ============================================================ real colvector _xpvc_BartlettCorr(real colvector stats_TR, real scalar dim_K, real scalar dim_p, real scalar dim_T, string scalar type) { real colvector corrected real scalar k, d, n_par, bart_factor corrected = stats_TR // Bartlett correction factor (Johansen 2002:Table 1) // Simplified: scale by T / (T - K*p - n_det) real scalar n_det if (type == "Case1") n_det = 0 else if (type == "Case2") n_det = 1 else if (type == "Case3") n_det = 1 else if (type == "Case4") n_det = 2 else if (type == "Case5") n_det = 2 else n_det = 1 n_par = dim_K * dim_p + n_det bart_factor = dim_T / (dim_T - n_par) // Apply Bartlett-type size correction // This is the Reinsel-Ahn correction (Reinsel & Ahn 1992) for (k = 1; k <= rows(stats_TR); k++) { corrected[k] = stats_TR[k] / bart_factor } return(corrected) } // ============================================================ // WALD JOINT TEST // Test linear restrictions R * vec(COEF) = r // ============================================================ struct _xpvc_WaldResult { real scalar chi2 real scalar df real scalar pval } struct _xpvc_WaldResult scalar _xpvc_WaldTest( real matrix COEF, real matrix VCOV, real matrix R, real colvector r) { struct _xpvc_WaldResult scalar res real colvector theta, Rtheta_r real matrix RVR theta = vec(COEF) Rtheta_r = R * theta - r RVR = R * VCOV * R' res.chi2 = Rtheta_r' * invsym(RVR) * Rtheta_r res.df = rows(R) res.pval = 1 - chi2(res.df, res.chi2) return(res) } // ============================================================ // VECM TO VAR: Deterministic terms conversion // Extract deterministic terms from VECM representation // ============================================================ // Compute the deterministic intercept/trend in VAR levels form // Given VECM: dY_t = alpha*beta'*Y_{t-1} + GAMMA*Z2_t + D*mu + u_t // Convert to levels VAR intercept nu = (I - A_1 - ... - A_p) * mu_0 void _xpvc_vec2var_det(real matrix alpha, real matrix beta, real matrix GAMMA_det, real matrix D_coef, real scalar dim_p, real scalar dim_K, real matrix nu_levels) { real matrix PI_K, I_K, A_sum real scalar j PI_K = alpha * beta[1..dim_K, .]' I_K = I(dim_K) // Compute sum of A matrices from vec2var A_sum = J(dim_K, dim_K, 0) real matrix A_full A_full = _xpvc_vec2var(alpha * beta', GAMMA_det, dim_p) for (j = 1; j <= dim_p; j++) { A_sum = A_sum + A_full[., (j-1)*dim_K+1..j*dim_K] } // nu = (I - sum(A_j)) * mu + PI * restricted_det if (cols(D_coef) > 0) { nu_levels = (I_K - A_sum) * D_coef } else { nu_levels = J(dim_K, 1, 0) } } // ============================================================ // STABILITY CHECK for VAR/VECM // Comprehensive stability assessment // ============================================================ struct _xpvc_StabilityResult { real colvector eigenmod // moduli of companion eigenvalues real scalar max_mod // maximum modulus real scalar is_stable // 1 if stable, 0 otherwise real scalar n_unit_roots // number of eigenvalues near unity } struct _xpvc_StabilityResult scalar _xpvc_stability( real matrix A, real scalar dim_p, | real scalar tol_unit) { struct _xpvc_StabilityResult scalar res real matrix C_comp real colvector eig_vals real matrix eig_vecs real scalar dim_K, j if (args() < 3) tol_unit = 0.001 dim_K = rows(A) // Build companion matrix C_comp = _xpvc_companion(A, dim_p) // Eigenvalue decomposition eigensystem(C_comp, eig_vecs, eig_vals) res.eigenmod = abs(eig_vals)' // transpose to column vector // Sort descending real colvector idx idx = order(-res.eigenmod, 1) res.eigenmod = res.eigenmod[idx] res.max_mod = max(res.eigenmod) res.is_stable = (res.max_mod < 1) // Count unit roots (eigenvalues within tol of 1.0) res.n_unit_roots = 0 for (j = 1; j <= rows(res.eigenmod); j++) { if (abs(res.eigenmod[j] - 1) < tol_unit) { res.n_unit_roots = res.n_unit_roots + 1 } } return(res) } // ============================================================ // LOG-LIKELIHOOD COMPUTATION for VECM // Johansen (1995:93, Eq.6.12) // ============================================================ real scalar _xpvc_loglik(real colvector lambda, real scalar dim_r, real scalar dim_T, real scalar dim_K) { real scalar ll, j // Concentrated log-likelihood // const + (-T/2) * sum_{j=1}^{r} ln(1 - lambda_j) // We just return the rank-dependent part ll = 0 for (j = 1; j <= dim_r; j++) { ll = ll + ln(1 - lambda[j]) } ll = -dim_T / 2 * ll return(ll) } // ============================================================ // CUMULATIVE SUM (CUSUM) RESIDUAL TEST // Brown, Durbin & Evans (1975) // ============================================================ void _xpvc_CUSUM(real matrix resid, real scalar dim_K, real matrix cusum_out, real matrix bounds_out) { real scalar dim_T, t real matrix cusum_t real colvector sigma_u dim_T = cols(resid) // Recursive residuals (simplified: just standardize by std dev) sigma_u = sqrt(diagonal(resid * resid' / dim_T)) // CUSUM: S_t = sum_{s=k+1}^{t} w_s / sigma cusum_out = J(dim_K, dim_T, 0) for (t = 2; t <= dim_T; t++) { cusum_out[., t] = cusum_out[., t-1] + resid[., t] :/ sigma_u } // 5% significance bounds: +/- a + 2*a*(t-k)/(T-k) where a = 0.948 real scalar a a = 0.948 bounds_out = J(2, dim_T, .) for (t = 1; t <= dim_T; t++) { bounds_out[1, t] = a * sqrt(dim_T) + 2 * a * (t - 1) / sqrt(dim_T) bounds_out[2, t] = -(a * sqrt(dim_T) + 2 * a * (t - 1) / sqrt(dim_T)) } } // ============================================================ // NORMALITY TEST (Doornik-Hansen / Jarque-Bera multivariate) // Doornik & Hansen (1994) // ============================================================ struct _xpvc_NormResult { real scalar chi2_jb // Jarque-Bera chi-squared real scalar df_jb // degrees of freedom real scalar pval_jb // p-value real colvector skewness // individual skewness real colvector kurtosis // individual kurtosis } struct _xpvc_NormResult scalar _xpvc_NormTest(real matrix resid, real scalar dim_K) { struct _xpvc_NormResult scalar res real scalar dim_T, k real matrix resid_t real colvector skew, kurt, jb_k real scalar chi2_total dim_T = cols(resid) resid_t = resid' // T x K skew = J(dim_K, 1, .) kurt = J(dim_K, 1, .) jb_k = J(dim_K, 1, .) for (k = 1; k <= dim_K; k++) { real colvector u real scalar m2, m3, m4, s, ke u = resid_t[., k] :- mean(resid_t[., k]) m2 = mean(u:^2) m3 = mean(u:^3) m4 = mean(u:^4) if (m2 > 0) { s = m3 / (m2^1.5) // skewness ke = m4 / (m2^2) - 3 // excess kurtosis } else { s = 0 ke = 0 } skew[k] = s kurt[k] = ke + 3 // raw kurtosis // Individual Jarque-Bera jb_k[k] = dim_T * (s^2 / 6 + ke^2 / 24) } chi2_total = sum(jb_k) res.chi2_jb = chi2_total res.df_jb = 2 * dim_K res.pval_jb = 1 - chi2(res.df_jb, chi2_total) res.skewness = skew res.kurtosis = kurt return(res) } // ============================================================ // LM RANK TEST — Breitung (2005:158) // Saikkonen (1999) / Breitung (2005) // ============================================================ struct _xpvc_LMrankResult { real colvector stats_LM real colvector pvals_LM real colvector r_H0 } struct _xpvc_LMrankResult scalar _xpvc_LMrank( real matrix R0, real matrix R1, pointer(real matrix) colvector L_alpha_oc, pointer(real matrix) colvector L_beta_oc, real matrix moments) { struct _xpvc_LMrankResult scalar res real scalar dim_T, dim_K, n_tests, r, idx real matrix U, W, UW, WWinv, UUinv, plus real scalar lambda_lm real colvector m, v dim_T = cols(R0) dim_K = rows(R0) n_tests = rows(L_alpha_oc) res.r_H0 = J(n_tests, 1, .) res.stats_LM = J(n_tests, 1, .) res.pvals_LM = J(n_tests, 1, .) // Restricted deterministic in auxiliary regression plus = J(0, dim_T, 0) if (rows(R1) > dim_K) { plus = R1[dim_K+1..rows(R1), .] } for (idx = 1; idx <= n_tests; idx++) { r = dim_K - cols(*L_alpha_oc[idx]) res.r_H0[idx] = r // Variables for auxiliary regression (Breitung 2005:158, Eq.13) U = (*L_alpha_oc[idx])' * R0 W = (*L_beta_oc[idx])' * R1[1..dim_K, .] if (rows(plus) > 0) { W = W \ plus } // LM statistic (Breitung 2005:158, Eq.13(II)) UW = U * W' WWinv = invsym(W * W') UUinv = invsym(U * U') lambda_lm = dim_T * trace(UW * WWinv * UW' * UUinv) res.stats_LM[idx] = lambda_lm } // Gamma-approximation p-values (Doornik 1998:578) m = moments[., 1] v = moments[., 2] for (idx = 1; idx <= n_tests; idx++) { if (m[idx] != . & v[idx] != . & v[idx] > 0) { real scalar shape, rate shape = m[idx]^2 / v[idx] rate = m[idx] / v[idx] res.pvals_LM[idx] = 1 - gammap(shape, rate * res.stats_LM[idx]) } else { res.pvals_LM[idx] = . } } return(res) } // ============================================================ // LR RESTRICTION TEST — Johansen (1995:126, Eq.8.17) // Test restrictions on alpha (weak exogeneity) // ============================================================ struct _xpvc_LRrestResult { real scalar stats real scalar pvals real scalar df } struct _xpvc_LRrestResult scalar _xpvc_LRrest( real colvector lambda, real colvector lambda_rest, real scalar dim_r, real scalar dim_T, real scalar dim_L) { struct _xpvc_LRrestResult scalar res real scalar j res.stats = 0 for (j = 1; j <= dim_r; j++) { res.stats = res.stats + ln((1 - lambda_rest[j]) / (1 - lambda[j])) } res.stats = dim_T * res.stats res.df = dim_r * dim_L res.pvals = 1 - chi2(res.df, res.stats) return(res) } // ============================================================ // ARCH-LM TEST — Luetkepohl (2005:576) // Multivariate ARCH test for residual heteroskedasticity // ============================================================ struct _xpvc_ARCHResult { real scalar LM_stat real scalar df real scalar pval real scalar lag_h } struct _xpvc_ARCHResult scalar _xpvc_ARCH( real matrix resid, real scalar dim_K, real scalar lag_h) { struct _xpvc_ARCHResult scalar res real scalar dim_T, n_vech, t, h, j real matrix vech_sel, u_sq, Y, Zh, Bh, eh, B0, e0 real matrix eh_cov, e0_cov real scalar TR, LM dim_T = cols(resid) // vech indices (lower triangular including diagonal) n_vech = dim_K * (dim_K + 1) / 2 // Squared residual vectors: vech(u_t u_t') u_sq = J(n_vech, dim_T, 0) for (t = 1; t <= dim_T; t++) { real colvector ut real matrix ut_outer ut = resid[., t] ut_outer = ut * ut' real scalar row_idx, ki, kj row_idx = 0 for (kj = 1; kj <= dim_K; kj++) { for (ki = kj; ki <= dim_K; ki++) { row_idx = row_idx + 1 u_sq[row_idx, t] = ut_outer[ki, kj] } } } // Auxiliary VAR(h) for u_sq Y = u_sq[., lag_h+1..dim_T] real scalar T_eff T_eff = cols(Y) // Stack lagged u_sq as regressors + constant Zh = J(1, T_eff, 1) // constant for (h = 1; h <= lag_h; h++) { Zh = Zh \ u_sq[., lag_h+1-h..dim_T-h] } // OLS VAR(h) Bh = Y * Zh' * invsym(Zh * Zh') eh = Y - Bh * Zh eh_cov = eh * eh' / T_eff // OLS VAR(0) — constant only real matrix Z0 Z0 = J(1, T_eff, 1) B0 = Y * Z0' * invsym(Z0 * Z0') e0 = Y - B0 * Z0 e0_cov = e0 * e0' / T_eff // LM statistic TR = trace(eh_cov * invsym(e0_cov)) LM = 0.5 * T_eff * n_vech - T_eff * TR res.LM_stat = LM res.lag_h = lag_h res.df = lag_h * n_vech^2 if (res.df > 0) { res.pval = 1 - chi2(res.df, LM) } else { res.pval = . } return(res) } // ============================================================ // SERIAL CORRELATION LM TEST — Luetkepohl (2005:173,347) // Breusch-Godfrey LM + Edgerton-Shukur F-test // ============================================================ struct _xpvc_SerialResult { real scalar LM_stat real scalar LM_pval real scalar F_stat real scalar F_pval real scalar lag_h } struct _xpvc_SerialResult scalar _xpvc_SerialCorr( real matrix resid, real matrix X_regs, real scalar dim_K, real scalar lag_h) { struct _xpvc_SerialResult scalar res real scalar dim_T, i, m, n_reg real matrix Y, Z, B, e, BR, eR real matrix e_cov, R_cov real scalar TR, LM, r_val, q_val, N_val, R_det, LMF dim_T = cols(resid) Y = resid Z = X_regs // Stack lagged residuals as additional regressors for (i = 1; i <= lag_h; i++) { real matrix u_lag u_lag = J(dim_K, i, 0), resid[., 1..dim_T-i] Z = Z \ u_lag } // Unrestricted OLS B = Y * Z' * invsym(Z * Z') e = Y - B * Z e_cov = e * e' / dim_T // Restricted OLS (original model only) BR = Y * X_regs' * invsym(X_regs * X_regs') eR = Y - BR * X_regs R_cov = eR * eR' / dim_T // Breusch-Godfrey LM TR = trace(invsym(R_cov) * e_cov) LM = dim_T * (dim_K - TR) res.LM_stat = LM res.lag_h = lag_h res.LM_pval = 1 - chi2(lag_h * dim_K^2, LM) // Edgerton-Shukur F-test m = dim_K * lag_h n_reg = rows(X_regs) if (dim_K^2 + m^2 - 5 > 0) { r_val = ((dim_K^2 * m^2 - 4) / (dim_K^2 + m^2 - 5))^0.5 } else { r_val = 1 } if (dim_K == 1) r_val = 1 q_val = 0.5 * dim_K * m - 1 N_val = dim_T - n_reg - m - 0.5 * (dim_K - m + 1) R_det = det(R_cov) / det(e_cov) if (R_det > 0 & R_det != .) { LMF = (R_det^(1/r_val) - 1) * (N_val * r_val - q_val) / (dim_K * m) res.F_stat = LMF res.F_pval = Ftail(lag_h * dim_K^2, N_val * r_val - q_val, LMF) } else { res.F_stat = . res.F_pval = . } return(res) } // ============================================================ // KILIAN (1998) BIAS-CORRECTION — Bootstrap-after-bootstrap // ============================================================ real matrix _xpvc_BaB(real matrix A_hat, real scalar dim_p, real matrix PSI) { real matrix A_bc, C_comp real colvector eig_vals real matrix eig_vecs real scalar m_hat, m_bc, delta_i, j real colvector deltas // Check stationarity of original estimate C_comp = _xpvc_companion(A_hat, dim_p) eigensystem(C_comp, eig_vecs, eig_vals) m_hat = max(abs(eig_vals)) if (m_hat >= 1) { // No bias-correction under non-stationarity return(A_hat) } // Geometric delta sequence: 1.0, 0.99, 0.98, ..., 0.0 deltas = (100::0) / 100 for (j = 2; j <= rows(deltas); j++) { deltas[j] = deltas[j-1] * deltas[j] } // Keep bias-corrected A within stationarity for (j = 1; j <= rows(deltas); j++) { delta_i = deltas[j] A_bc = A_hat - delta_i * PSI C_comp = _xpvc_companion(A_bc, dim_p) eigensystem(C_comp, eig_vecs, eig_vals) m_bc = max(abs(eig_vals)) if (m_bc < 1) return(A_bc) } return(A_hat) } // ============================================================ // CONDITIONAL → FULL VECM TRANSFORM // Johansen (1995:122, Th.8.3) // ============================================================ void _xpvc_con2vec( real matrix alpha_y, real matrix GAMMA_c, real matrix resid_c, real matrix OMEGA_c, real matrix SIGMA_c, real matrix GAMMA_x, real matrix resid_x, real matrix OMEGA_xx, real matrix SIGMA_xx, real matrix LAMBDA, real matrix beta, real scalar dim_K, real scalar dim_L, real scalar dim_r, real matrix alpha_out, real matrix PI_out, real matrix GAMMA_out, real matrix OMEGA_out, real matrix SIGMA_out, real matrix resid_out) { real matrix GAMMA_y_tr, resid_y real matrix alpha_x, OMEGA_yx, SIGMA_yx real scalar dim_T dim_T = min((cols(resid_c), cols(resid_x))) // Transform short-run (Johansen 1995:122, Th.8.3) GAMMA_y_tr = GAMMA_c + LAMBDA * GAMMA_x resid_y = resid_c[., 1..dim_T] + LAMBDA * resid_x[., 1..dim_T] // Alpha under weak exogeneity alpha_x = J(dim_L, dim_r, 0) // Covariance transformations OMEGA_yx = LAMBDA * OMEGA_xx SIGMA_yx = LAMBDA * SIGMA_xx // Stack to full system alpha_out = alpha_y \ alpha_x PI_out = alpha_out * beta' GAMMA_out = GAMMA_y_tr \ GAMMA_x resid_out = resid_y \ resid_x[., 1..dim_T] OMEGA_out = (OMEGA_c + LAMBDA * OMEGA_yx', OMEGA_yx) \ /// (OMEGA_yx', OMEGA_xx) SIGMA_out = (SIGMA_c + LAMBDA * SIGMA_yx', SIGMA_yx) \ /// (SIGMA_yx', SIGMA_xx) } // ============================================================ // CONDITIONAL → FULL VAR TRANSFORM // Luetkepohl (2005:392) // ============================================================ void _xpvc_con2var( real matrix A_c, real matrix resid_c, real matrix OMEGA_c, real matrix SIGMA_c, real matrix A_xx, real matrix resid_x, real matrix OMEGA_xx, real matrix SIGMA_xx, real matrix LAMBDA, real scalar dim_K, real scalar dim_L, real matrix A_out, real matrix OMEGA_out, real matrix SIGMA_out, real matrix resid_out) { real matrix A_y_tr, resid_y real matrix OMEGA_yx, SIGMA_yx real scalar dim_T dim_T = min((cols(resid_c), cols(resid_x))) // Transform coefficient matrix A_y_tr = A_c + LAMBDA * A_xx resid_y = resid_c[., 1..dim_T] + LAMBDA * resid_x[., 1..dim_T] // Covariance OMEGA_yx = LAMBDA * OMEGA_xx SIGMA_yx = LAMBDA * SIGMA_xx // Stack A_out = A_y_tr \ A_xx resid_out = resid_y \ resid_x[., 1..dim_T] OMEGA_out = (OMEGA_c + LAMBDA * OMEGA_yx', OMEGA_yx) \ /// (OMEGA_yx', OMEGA_xx) SIGMA_out = (SIGMA_c + LAMBDA * SIGMA_yx', SIGMA_yx) \ /// (SIGMA_yx', SIGMA_xx) } // ============================================================ // VECM → VMA: Granger Representation Theorem // Luetkepohl (2005:252, Eq.6.3.12) // ============================================================ struct _xpvc_VMAResult { real matrix XI // long-run multiplier (K x K) real matrix UPSILON // structural long-run multiplier (K x S) } struct _xpvc_VMAResult scalar _xpvc_vec2vma( real matrix GAMMA, real matrix alpha_oc, real matrix beta_oc, real scalar dim_p, real matrix B) { struct _xpvc_VMAResult scalar res real scalar dim_K, j real matrix GAMMA_1, A_sum dim_K = rows(GAMMA) // Cumulative GAMMA(1) = I - sum(GAMMA_j) A_sum = J(dim_K, dim_K, 0) if (dim_p > 1) { real scalar n_gamma_cols n_gamma_cols = dim_K * (dim_p - 1) if (n_gamma_cols <= cols(GAMMA)) { for (j = 1; j <= dim_p - 1; j++) { A_sum = A_sum + GAMMA[., (j-1)*dim_K+1..j*dim_K] } } } GAMMA_1 = I(dim_K) - A_sum // Long-run multiplier matrix (Luetkepohl 2005:252, Eq.6.3.12) res.XI = beta_oc * invsym(alpha_oc' * GAMMA_1 * beta_oc) * alpha_oc' // Structural long-run multiplier res.UPSILON = res.XI * B return(res) } // ============================================================ // PERSISTENCE PROFILES — System-wide shock // Pesaran, Shin (1996:125, Eq.11/13) // ============================================================ void _xpvc_PP_system( real matrix A, real scalar dim_p, real matrix beta, real matrix SIGMA, real scalar dim_r, real scalar n_ahead, real matrix hz_out) { real scalar dim_K, n, j real matrix C_comp, PHI_n, B_n, Hz_0, G real colvector g dim_K = rows(A) C_comp = _xpvc_companion(A, dim_p) // Compute persistence profiles hz_out = J(dim_r * dim_r, n_ahead + 1, .) PHI_n = I(dim_K * dim_p) real matrix Hz_init Hz_init = J(dim_r, dim_r, 0) for (n = 0; n <= n_ahead; n++) { B_n = PHI_n[1..dim_K, 1..dim_K] // Unscaled PP (Pesaran,Shin 1996:125, Eq.11) real matrix Hz_n Hz_n = beta' * B_n * SIGMA * B_n' * beta // Scale at impact (n=0) if (n == 0) Hz_init = Hz_n // Scaling matrix from impact coefficients g = J(dim_r, 1, 1) for (j = 1; j <= dim_r; j++) { if (Hz_init[j, j] > 0) g[j] = 1 / sqrt(Hz_init[j, j]) } G = diag(g) // Scaled PP (Pesaran,Shin 1996:125, Eq.13) real matrix hz_n hz_n = G * Hz_n * G hz_out[., n+1] = vec(hz_n) PHI_n = PHI_n * C_comp } } // ============================================================ // PERSISTENCE PROFILES — Variable-specific shock // Pesaran, Shin (1996:122, Eq.8) // ============================================================ void _xpvc_PP_variable( real matrix A, real scalar dim_p, real matrix beta, real matrix U_P, real matrix shock, real scalar dim_r, real scalar n_ahead, real matrix psi_out) { real scalar dim_K, n real matrix C_comp, PHI_n, B_n dim_K = rows(A) C_comp = _xpvc_companion(A, dim_p) // PP: psi_out is (dim_r * n_shocks) x (n_ahead+1) real scalar n_shocks n_shocks = cols(shock) psi_out = J(dim_r * n_shocks, n_ahead + 1, .) PHI_n = I(dim_K * dim_p) for (n = 0; n <= n_ahead; n++) { B_n = PHI_n[1..dim_K, 1..dim_K] // PP (Pesaran,Shin 1996:122, Eq.8) real matrix PSI_n PSI_n = beta' * B_n * U_P * shock psi_out[., n+1] = vec(PSI_n) PHI_n = PHI_n * C_comp } } // ============================================================ // PROXY SVAR IDENTIFICATION — Mertens-Ravn (2013) // Jentsch-Lunsford (2019) // ============================================================ struct _xpvc_ProxySVARResult { real matrix B1 // structural impact matrix (K x S) real matrix shocks // estimated structural shocks (S x T) real matrix F_stats // F-test for weak instruments } struct _xpvc_ProxySVARResult scalar _xpvc_idIV_MR( real matrix u, real matrix m, real matrix SIGMA_uu) { struct _xpvc_ProxySVARResult scalar res real scalar dim_T, dim_K, dim_S real matrix SIGMA_mm, SIGMA_um, SIGMA_mu real matrix B21B11inv, M1_tmp, bigZ, M2_tmp real matrix B12B12, B11B11, B22B22, B12B22inv real matrix M3_tmp, S1S1, B11, B1 dim_T = cols(u) dim_K = rows(u) dim_S = rows(m) SIGMA_mm = m * m' / dim_T SIGMA_um = u * m' / dim_T SIGMA_mu = SIGMA_um' // Mertens-Ravn (2013:1244) real matrix SIGMA_mu11, SIGMA_mu12 real matrix SIGMA_uu11, SIGMA_uu21, SIGMA_uu22 SIGMA_mu11 = SIGMA_mu[1..dim_S, 1..dim_S] SIGMA_mu12 = SIGMA_mu[1..dim_S, dim_S+1..dim_K] SIGMA_uu11 = SIGMA_uu[1..dim_S, 1..dim_S] SIGMA_uu21 = SIGMA_uu[dim_S+1..dim_K, 1..dim_S] SIGMA_uu22 = SIGMA_uu[dim_S+1..dim_K, dim_S+1..dim_K] B21B11inv = (invsym(SIGMA_mu11) * SIGMA_mu12)' // Eq.10 M1_tmp = B21B11inv * SIGMA_uu21' bigZ = SIGMA_uu22 - M1_tmp - M1_tmp' + B21B11inv * SIGMA_uu11 * B21B11inv' M2_tmp = SIGMA_uu21 - B21B11inv * SIGMA_uu11 B12B12 = M2_tmp' * invsym(bigZ) * M2_tmp B11B11 = SIGMA_uu11 - B12B12 B22B22 = SIGMA_uu22 - B21B11inv * B11B11 * B21B11inv' B12B22inv = (SIGMA_uu21' - B11B11 * B21B11inv') * invsym(B22B22) M3_tmp = I(dim_S) - B12B22inv * B21B11inv S1S1 = M3_tmp * B11B11 * M3_tmp' // Eq.13 B11 = invsym(M3_tmp) * cholesky(S1S1)' // Eq.14 (lower triangular) res.B1 = B11 \ (B21B11inv * B11) // Eq.9+10 // Shocks and F-stats real matrix Pi, PSIinv Pi = invsym(SIGMA_uu) * SIGMA_um PSIinv = invsym(cholesky(SIGMA_mu * Pi)) res.shocks = (Pi * PSIinv)' * u // F-test for weak instrument (Lunsford 2016:15) real matrix SIGMA_ee SIGMA_ee = (m - Pi' * u) * (m - Pi' * u)' res.F_stats = ((dim_T - dim_K) / dim_K) * (SIGMA_mm - SIGMA_ee) :/ SIGMA_ee return(res) } // Jentsch-Lunsford (2019:34) proxy SVAR struct _xpvc_ProxySVARResult scalar _xpvc_idIV_JL( real matrix u, real matrix m, real matrix SIGMA_uu) { struct _xpvc_ProxySVARResult scalar res real scalar dim_T, dim_K, dim_S real matrix SIGMA_mm, SIGMA_um, SIGMA_mu real matrix Pi, PSIinv dim_T = cols(u) dim_K = rows(u) dim_S = rows(m) SIGMA_mm = m * m' / dim_T SIGMA_um = u * m' / dim_T SIGMA_mu = SIGMA_um' Pi = invsym(SIGMA_uu) * SIGMA_um PSIinv = invsym(cholesky(SIGMA_mu * Pi)) res.B1 = SIGMA_um * PSIinv res.shocks = (Pi * PSIinv)' * u // F-test real matrix SIGMA_ee SIGMA_ee = (m - Pi' * u) * (m - Pi' * u)' res.F_stats = ((dim_T - dim_K) / dim_K) * (SIGMA_mm - SIGMA_ee) :/ SIGMA_ee return(res) } // ============================================================ // SCORING ALGORITHM for SVEC/GRT identification // Amisano,Giannini (1997:57) / Luetkepohl,Kraetzig (2004:173) // ============================================================ struct _xpvc_SVECResult { real matrix SR // short-run impact matrix real matrix LR // long-run impact matrix real scalar iter // iterations to convergence } struct _xpvc_SVECResult scalar _xpvc_idGRT_scoring( real matrix OMEGA, real matrix XI, real scalar dim_T, real matrix LR_mask, real matrix SR_mask, real scalar max_iter, real scalar conv_crit, real scalar maxls) { struct _xpvc_SVECResult scalar res real scalar dim_K, K2, l, iters real matrix IK, IK2, Kkk real matrix R_B, R_C1, R_mat, Sb, S_mat real colvector s_vec, gamma, vecab real matrix A_mat, B_mat, Binv, Btinv, BinvA real matrix infvecab, infgamma, infgammainv real colvector scorevecBinvA, scorevecAB, scoregamma, direction real scalar lambda_step, cvcrit, length_dir real colvector z dim_K = rows(OMEGA) K2 = dim_K^2 IK = I(dim_K) IK2 = I(K2) // Commutation matrix Kkk = J(K2, K2, 0) real scalar ci, ri for (ci = 1; ci <= dim_K; ci++) { for (ri = 1; ri <= dim_K; ri++) { Kkk[(ci-1)*dim_K + ri, (ri-1)*dim_K + ci] = 1 } } // Build restriction matrices from mask // SR_mask: 0 = restricted to zero, . = free // LR_mask: 0 = restricted to zero, . = free real scalar n_free_B, n_free_C, idx_r R_B = J(0, K2, 0) for (idx_r = 1; idx_r <= K2; idx_r++) { if (SR_mask[mod(idx_r-1, dim_K)+1, floor((idx_r-1)/dim_K)+1] == 0) { real rowvector row_sel row_sel = J(1, K2, 0) row_sel[idx_r] = 1 R_B = R_B \ row_sel } } R_C1 = J(0, K2, 0) for (idx_r = 1; idx_r <= K2; idx_r++) { if (LR_mask[mod(idx_r-1, dim_K)+1, floor((idx_r-1)/dim_K)+1] == 0) { real rowvector row_sel2 row_sel2 = J(1, K2, 0) row_sel2[idx_r] = 1 R_C1 = R_C1 \ row_sel2 } } R_C1 = R_C1 * (XI # IK) if (rows(R_B) == 0) R_mat = R_C1 else R_mat = R_C1 \ R_B // S matrix (explicit form) via orthogonal complement Sb = _xpvc_oc(R_mat') S_mat = J(K2, cols(Sb), 0) \ Sb l = cols(S_mat) s_vec = vec(IK) \ J(K2, 1, 0) // Initialize gamma = rnormal(l, 1, 0, 1) // Scoring algorithm iters = 0 cvcrit = conv_crit + 1 while (cvcrit > conv_crit) { z = gamma vecab = S_mat * gamma + s_vec A_mat = rowshape(vecab[1..K2], dim_K) B_mat = rowshape(vecab[K2+1..2*K2], dim_K) Binv = invsym(B_mat) Btinv = invsym(B_mat') BinvA = Binv * A_mat // Information matrix real matrix inf1, inf2, inf3 inf1 = (invsym(BinvA) # Btinv) \ (-IK # Btinv) inf2 = IK2 + Kkk inf3 = (BinvA' # Binv, -IK # Binv) infvecab = dim_T * (inf1 * inf2 * inf3) infgamma = S_mat' * infvecab * S_mat infgammainv = invsym(infgamma) // Score scorevecBinvA = dim_T * vec(invsym(BinvA')) - dim_T * ((OMEGA # IK) * vec(BinvA)) real matrix scoreAB_mat scoreAB_mat = (IK # Btinv) \ (-(BinvA # Btinv)) scorevecAB = scoreAB_mat * scorevecBinvA scoregamma = S_mat' * scorevecAB direction = infgammainv * scoregamma length_dir = max(abs(direction)) if (length_dir > maxls) lambda_step = maxls / length_dir else lambda_step = 1 gamma = gamma + lambda_step * direction iters = iters + 1 z = z - gamma cvcrit = max(abs(z)) if (iters >= max_iter) break } // Final B matrix vecab = S_mat * gamma + s_vec A_mat = rowshape(vecab[1..K2], dim_K) B_mat = rowshape(vecab[K2+1..2*K2], dim_K) // Normalize sign real matrix BinvA_final BinvA_final = invsym(A_mat) * B_mat real scalar k for (k = 1; k <= dim_K; k++) { if (BinvA_final[k, k] < 0) B_mat[., k] = -B_mat[., k] } res.SR = B_mat res.LR = XI * B_mat res.iter = iters - 1 // Clean near-zero elements for (ri = 1; ri <= dim_K; ri++) { for (ci = 1; ci <= dim_K; ci++) { if (abs(res.LR[ri, ci]) < 1e-8) res.LR[ri, ci] = 0 if (abs(res.SR[ri, ci]) < 1e-8) res.SR[ri, ci] = 0 } } return(res) } // ============================================================ // BAI-NG (2002) PANEL INFORMATION CRITERIA — Full version // PC, IC, IPC variants // ============================================================ real matrix _xpvc_PIC(real scalar dim_KN, real scalar dim_T, real scalar k, real scalar Vk0, real scalar Vkmax0) { real scalar n_vals, C2_NT, p1, p2, p3, alphaT, ip3 real matrix result n_vals = dim_KN * dim_T C2_NT = min((dim_KN, dim_T)) // Bai-Ng (2002:201, Eq.9) p1 = (dim_KN + dim_T) / n_vals * ln(n_vals / (dim_KN + dim_T)) p2 = (dim_KN + dim_T) / n_vals * ln(C2_NT) p3 = ln(C2_NT) / C2_NT real rowvector penalties penalties = (p1, p2, p3) // PC criteria real rowvector PC, IC PC = J(1, 3, Vk0) + k * Vkmax0 * penalties IC = J(1, 3, ln(Vk0)) + k * penalties // IPC criteria (Bai 2004:145, Eq.12) alphaT = dim_T / (4 * ln(ln(dim_T))) ip3 = (dim_KN + dim_T - k) / n_vals * ln(n_vals) real rowvector IPC IPC = J(1, 3, Vk0) + k * Vkmax0 * alphaT * (p1, p2, ip3) result = PC \ IC \ IPC return(result) } // ============================================================ // COEFFICIENT ACCUMULATION A(1) — Luetkepohl (2005:289) // ============================================================ real matrix _xpvc_accum(real matrix A, real scalar dim_p) { real scalar dim_K, j real matrix result dim_K = rows(A) result = I(dim_K) for (j = 1; j <= dim_p; j++) { result = result - A[., (j-1)*dim_K+1..j*dim_K] } return(result) } // ============================================================ // GLS DETRENDING with Q-TRANSFORMATION — Option 2 // Saikkonen, Luetkepohl (2000:439, Eq.3.4) // ============================================================ struct _xpvc_GLSresult scalar _xpvc_GLStrend_Q( real matrix y, real scalar dim_p, real matrix OMEGA, real matrix A_slope, real matrix D, real matrix alpha, real matrix alpha_oc) { struct _xpvc_GLSresult scalar res real scalar dim_T, dim_K, dim_n, dim_r, j real matrix OMEGAinv, Aj, Yj, Yz, Q dim_T = cols(y) dim_K = rows(y) dim_n = rows(D) dim_r = cols(alpha) OMEGAinv = invsym(OMEGA) // Build Q transformation matrix (Saikkonen,Luetkepohl 2000:438, Eq.3.2) real matrix Q1_left, Q2_right Q1_left = J(dim_K, 0, 0) Q2_right = J(dim_K, 0, 0) if (dim_r > 0) { real matrix Q1_tmp Q1_tmp = alpha' * OMEGAinv * alpha // Eigendecomposition for matrix square root inverse real matrix Q1_vecs real colvector Q1_vals eigensystem(Q1_tmp, Q1_vecs, Q1_vals) real matrix Q1_eival Q1_eival = diag(1 :/ sqrt(Re(Q1_vals))) Q1_left = OMEGAinv * alpha * /// Re(Q1_vecs) * Q1_eival * Re(Q1_vecs)' } if (dim_K - dim_r > 0) { real matrix Q2_tmp Q2_tmp = alpha_oc' * OMEGA * alpha_oc real matrix Q2_vecs real colvector Q2_vals eigensystem(Q2_tmp, Q2_vecs, Q2_vals) real matrix Q2_eival Q2_eival = diag(1 :/ sqrt(Re(Q2_vals))) Q2_right = alpha_oc * /// Re(Q2_vecs) * Q2_eival * Re(Q2_vecs)' } Q = Q1_left, Q2_right // Quasi-difference the data Yz = y for (j = 1; j <= dim_p; j++) { Aj = A_slope[., (j-1)*dim_K+1..j*dim_K] Yj = J(dim_K, j, 0), y[., 1..dim_T-j] Yz = Yz - Aj * Yj } // OLS on Q-transformed model (Saikkonen,Luetkepohl 2000:439, Eq.3.4) // Build QD and QY per observation, accumulate real matrix XtX, XtY real colvector mu_vec XtX = J(dim_K * dim_n, dim_K * dim_n, 0) XtY = J(dim_K * dim_n, 1, 0) real scalar t, ki, kj for (t = 1; t <= dim_T; t++) { real colvector d_t, qyz_t real matrix qd, dd d_t = D[., t] qyz_t = Q' * Yz[., t] dd = d_t * d_t' for (ki = 1; ki <= dim_K; ki++) { for (kj = 1; kj <= dim_K; kj++) { XtX[(ki-1)*dim_n+1..ki*dim_n, (kj-1)*dim_n+1..kj*dim_n] = /// XtX[(ki-1)*dim_n+1..ki*dim_n, (kj-1)*dim_n+1..kj*dim_n] + /// dd * (Q' * I(dim_K))[ki, kj] // simplified } XtY[(ki-1)*dim_n+1..ki*dim_n] = /// XtY[(ki-1)*dim_n+1..ki*dim_n] + d_t * qyz_t[ki] } } mu_vec = invsym(XtX) * XtY res.MU = rowshape(mu_vec, dim_K) res.x = y - res.MU * D return(res) } end