*! _dnqr_init.mata version 1.0.0 27may2026 *! Mata function library for the dnqr package. Loaded once per Stata *! session by the _dnqr_mata loader; do not call directly. *! *! Author : Dr Merwan Roudane version 13.0 mata: mata set matastrict off // --------------------------------------------------------------------- real matrix dnqr_rowstd(real matrix W) { real colvector rs real matrix Wo Wo = W rs = rowsum(Wo) rs = rs + (rs:==0) return(Wo :/ rs) } // --------------------------------------------------------------------- real matrix dnqr_panel_to_NT(real colvector v, real colvector id, real colvector t, real scalar N, real scalar T) { real matrix M real scalar i M = J(N, T, .) for (i=1; i<=rows(v); i++) M[id[i], t[i]] = v[i] return(M) } // --------------------------------------------------------------------- real matrix dnqr_make_F(real matrix Fmat, real scalar N, real scalar plag) { real scalar T, m, k, j real matrix Fdesign, Fblock T = rows(Fmat) m = cols(Fmat) if (plag<0) plag = 0 if (T-plag <= 0) _error(3300, "common-factor lag exceeds T") Fblock = J(T-plag, 0, .) for (k=0; k<=plag; k++) { Fblock = Fblock, Fmat[(plag+1-k)..(T-k), .] } Fdesign = J(N*(T-plag), cols(Fblock), .) for (j=1; j<=T-plag; j++) { Fdesign[((j-1)*N+1)::(j*N), .] = J(N,1,1) * Fblock[j, .] } return(Fdesign) } // --------------------------------------------------------------------- real matrix dnqr_make_Z(real matrix Zmat, real scalar Teff) { real matrix Zdesign real scalar j, N N = rows(Zmat) Zdesign = J(N*Teff, cols(Zmat), .) for (j=1; j<=Teff; j++) { Zdesign[((j-1)*N+1)::(j*N), .] = Zmat } return(Zdesign) } // --------------------------------------------------------------------- real scalar dnqr_powell_band(real scalar tau, real scalar n, string scalar btype, real scalar scale) { real scalar hf if (btype == "HS") { hf = n^(-1/3) * (invnormal(0.975))^(2/3) * (1.5 * (normalden(invnormal(tau)))^2 / (2 * (invnormal(tau))^2 + 1))^(1/3) } else { hf = n^(-1/5) * (4.5 * (normalden(invnormal(tau)))^4 / ((2*(invnormal(tau))^2 + 1)^2))^(1/5) } hf = scale * hf if (tau+hf>=1 | tau-hf<=0) hf = min((abs(tau-0.0001),abs(1-tau-0.0001))) return(hf) } // --------------------------------------------------------------------- real rowvector dnqr_qreg_via_st(real colvector y, real matrix X, real scalar tau) { real rowvector b real scalar k, j, rc, dummy string rowvector xn stata("preserve") stata("quietly drop _all") dummy = st_addobs(rows(y)) dummy = st_addvar("double", "__dnqr_y") st_store(., "__dnqr_y", y) k = cols(X) xn = J(1, k, "") for (j=1; j<=k; j++) { xn[j] = sprintf("__dnqr_x%g", j) dummy = st_addvar("double", xn[j]) st_store(., xn[j], X[., j]) } rc = _stata("_dnqr_silent_qreg __dnqr_y " + invtokens(xn) + ", quantile(" + strofreal(tau) + ")") if (rc != 0 | st_numscalar("c(rc)") != 0) { stata("restore") _error(3498, "qreg failed inside IVQR grid") } b = st_matrix("e(b)") stata("restore") return(b) } // --------------------------------------------------------------------- real rowvector dnqr_ivqr_grid(real colvector y, real colvector D, real matrix X, real matrix Mz, real colvector grid, real scalar tau) { real scalar i, L, gmin, alpha real rowvector b real colvector gnorms L = cols(Mz) gnorms = J(rows(grid), 1, .) for (i=1; i<=rows(grid); i++) { alpha = grid[i] b = dnqr_qreg_via_st(y :- alpha:*D, (Mz, X), tau) gnorms[i] = sqrt(sum(b[1..L]:^2)) } gmin = min(gnorms) for (i=1; i<=rows(grid); i++) { if (gnorms[i]==gmin) return((grid[i], gmin)) } return((grid[1], gnorms[1])) } // --------------------------------------------------------------------- real matrix dnqr_powell_vcov(real colvector y, real colvector D, real matrix X, real matrix Mz, real colvector resid, real scalar tau, string scalar btype, real scalar scale) { real scalar n, hf, den, qhi, qlo real matrix XD, XZ, J1, S, V n = rows(y) XD = (J(n,1,1), D, X) XZ = (J(n,1,1), Mz, X) hf = dnqr_powell_band(tau, n, btype, scale) qhi = dnqr_quantile(resid, tau+hf) qlo = dnqr_quantile(resid, tau-hf) if (qhi - qlo <= 0) den = 1e-6 else den = 2*hf / (qhi - qlo) J1 = (XD' * XZ) * den S = cross(XZ, XZ) V = invsym(J1 * invsym(S) * J1') * tau * (1-tau) return(V) } // --------------------------------------------------------------------- real scalar dnqr_quantile(real colvector x, real scalar p) { real colvector xs real scalar n, idx xs = sort(x, 1) n = rows(xs) idx = p * (n + 1) idx = max((1, min((n, idx)))) return(xs[round(idx)]) } // --------------------------------------------------------------------- real scalar dnqr_goodfit(real colvector r_full, real colvector r_null, real scalar tau) { real scalar V1, V0 V1 = sum(r_full :* (tau :- (r_full :< 0))) V0 = sum(r_null :* (tau :- (r_null :< 0))) if (V0==0) return(.) return(1 - V1/V0) } // --------------------------------------------------------------------- real matrix dnqr_impulse(real matrix G, real colvector v0, real scalar H) { real matrix IRF real scalar h IRF = J(rows(G), H+1, .) IRF[., 1] = v0 for (h=2; h<=H+1; h++) IRF[., h] = G * IRF[., h-1] return(IRF) } // ===================================================================== // Mata routines specific to the nqar / dnqr engines (recoding, design) // ===================================================================== real matrix _nqar_recode(real colvector v, real colvector u) { real scalar i, n real matrix idx real colvector r n = rows(v) r = J(n, 1, .) for (i=1; i<=n; i++) { idx = selectindex(u:==v[i]) r[i] = idx[1] } return(r) } real matrix _nqar_panel_to_NT(real colvector v, real colvector id2, real colvector t2, real scalar N, real scalar T) { real matrix M real scalar i M = J(N, T, .) for (i=1; i<=rows(v); i++) M[id2[i], t2[i]] = v[i] return(M) } real matrix _nqar_collapse_Z(real matrix Zlong, real colvector id2, real scalar N) { real matrix Z real scalar i Z = J(N, cols(Zlong), .) for (i=1; i<=rows(Zlong); i++) { if (Z[id2[i], 1] == .) Z[id2[i], .] = Zlong[i, .] } return(Z) } real matrix _nqar_collapse_F(real matrix Flong, real colvector t2, real scalar T) { real matrix F real scalar i F = J(T, cols(Flong), .) for (i=1; i<=rows(Flong); i++) { if (F[t2[i], 1] == .) F[t2[i], .] = Flong[i, .] } return(F) } real matrix _nqar_build_F(real matrix Fsub, real scalar N, real scalar plag) { real scalar Teff, m, k, j real matrix block, design Teff = rows(Fsub) - plag m = cols(Fsub) block = J(Teff, 0, .) for (k=0; k<=plag; k++) { block = block, Fsub[(plag+1-k)..(Teff+plag-k), .] } design = J(N*Teff, cols(block), .) for (j=1; j<=Teff; j++) { design[((j-1)*N+1)::(j*N), .] = J(N,1,1) * block[j, .] } return(design) } real matrix _nqar_expand_Z(real matrix Zmat, real scalar Teff) { real matrix design real scalar j, N N = rows(Zmat) if (cols(Zmat)==0) return(J(N*Teff, 0, .)) design = J(N*Teff, cols(Zmat), .) for (j=1; j<=Teff; j++) { design[((j-1)*N+1)::(j*N), .] = Zmat } return(design) } // Powell SE for the no-IV case real matrix _nqar_powell_vcov(real colvector y, real matrix X, real colvector r, real scalar tau, string scalar btype, real scalar bs) { real scalar n, hf, den, qhi, qlo real matrix Xc, V, A n = rows(y) Xc = (J(n,1,1), X) hf = dnqr_powell_band(tau, n, btype, bs) qhi = dnqr_quantile(r, tau+hf) qlo = dnqr_quantile(r, tau-hf) if (qhi - qlo <= 0) den = 1e-6 else den = 2*hf / (qhi - qlo) A = cross(Xc, Xc) :* den V = invsym(A * invsym(cross(Xc, Xc)) * A) :* (tau*(1-tau)) return(V) } // ===================================================================== // _nqar_engine : main NQAR estimator // ===================================================================== void _nqar_engine(string scalar depv, string scalar zlist, string scalar flist, string scalar Wname, string scalar pvar, string scalar tvar, string scalar touse, real scalar plag, string scalar qstr, string scalar btype, real scalar bscale, real scalar rowstd, real scalar lev, string scalar Bname, string scalar SEname, string scalar Tname, string scalar Pname, string scalar Loname, string scalar Hiname, string scalar Qname, string scalar Dname, string scalar Bfirst, string scalar Vfirst) { real matrix W, Ymat, Zmat, Fmat, X, Zd, Fd real matrix Ymatl, WYmatl, Zlong, Flong, Fsub real matrix B, SE, Tv, Pv, Lo, Hi real matrix V real colvector yv, idv, tv, id2, t2, ylag, wylag, taus real colvector uid, ut, bvec, resid real rowvector b real scalar N, T, Teff, q, nq, k, j, tstart, nb, netd, zlev, kreg string rowvector rnames, zn, fn string matrix rs, cs W = st_matrix(Wname) N = rows(W) if (cols(W) != N) _error(3200, "W must be square N x N") if (rowstd) W = dnqr_rowstd(W) idv = st_data(., pvar, touse) tv = st_data(., tvar, touse) yv = st_data(., depv, touse) uid = uniqrows(idv) ut = uniqrows(tv) if (rows(uid) != N) { printf("{err}rows(W) = %g but panel restricts to %g units\n", N, rows(uid)) _error(3200, "network dimension mismatch") } T = rows(ut) id2 = _nqar_recode(idv, uid) t2 = _nqar_recode(tv, ut) Ymat = _nqar_panel_to_NT(yv, id2, t2, N, T) if (zlist != "") { Zlong = st_data(., zlist, touse) Zmat = _nqar_collapse_Z(Zlong, id2, N) } else Zmat = J(N, 0, .) if (flist != "") { Flong = st_data(., flist, touse) Fmat = _nqar_collapse_F(Flong, t2, T) } else Fmat = J(T, 0, .) tstart = max((2, plag+2)) Teff = T - tstart + 1 if (Teff < 2) _error(3300, "not enough time periods after lags") Ymatl = Ymat[., (tstart-1)..(T-1)] WYmatl = W * Ymatl yv = vec(Ymat[., tstart..T]) ylag = vec(Ymatl) wylag = vec(WYmatl) Zd = _nqar_expand_Z(Zmat, Teff) if (cols(Fmat) > 0) { Fsub = Fmat[(tstart-plag)..T, .] Fd = _nqar_build_F(Fsub, N, plag) } else Fd = J(rows(yv), 0, .) X = wylag, ylag, Zd, Fd kreg = cols(X) rnames = ("_cons", "WY_L1", "Y_L1") if (cols(Zmat) > 0) rnames = rnames, tokens(zlist) if (cols(Fmat) > 0) { fn = tokens(flist) for (k=0; k<=plag; k++) { for (j=1; j<=cols(Fmat); j++) { if (k==0) rnames = rnames, fn[j] else rnames = rnames, sprintf("%s_L%g", fn[j], k) } } } nb = cols(rnames) taus = strtoreal(tokens(qstr))' nq = rows(taus) B = J(nb, nq, .); SE = J(nb, nq, .) Tv = J(nb, nq, .); Pv = J(nb, nq, .) Lo = J(nb, nq, .); Hi = J(nb, nq, .) zlev = invnormal(1 - (1-lev/100)/2) real matrix V_first, b_first V_first = J(nb, nb, .) b_first = J(1, nb, .) for (q=1; q<=nq; q++) { b = dnqr_qreg_via_st(yv, X, taus[q]) bvec = (b[kreg+1] \ b[1..kreg]') resid = yv :- (J(rows(yv),1,1), X) * bvec V = _nqar_powell_vcov(yv, X, resid, taus[q], btype, bscale) B[., q] = bvec SE[., q] = sqrt(diagonal(V)) Tv[., q] = B[., q] :/ SE[., q] Pv[., q] = 2 :* (1 :- normal(abs(Tv[., q]))) Lo[., q] = B[., q] :- zlev :* SE[., q] Hi[., q] = B[., q] :+ zlev :* SE[., q] if (q==1) { // store first-quantile b and V for ereturn-post b_first = bvec' V_first = V } } netd = sum(W:>0) / (N*N - N) st_numscalar(Dname, netd) rs = J(nb, 1, ""), rnames' cs = J(nq, 1, ""), strofreal(taus, "%6.3f") st_matrix(Bname, B); st_matrixrowstripe(Bname, rs); st_matrixcolstripe(Bname, cs) st_matrix(SEname, SE); st_matrixrowstripe(SEname, rs); st_matrixcolstripe(SEname, cs) st_matrix(Tname, Tv); st_matrixrowstripe(Tname, rs); st_matrixcolstripe(Tname, cs) st_matrix(Pname, Pv); st_matrixrowstripe(Pname, rs); st_matrixcolstripe(Pname, cs) st_matrix(Loname, Lo); st_matrixrowstripe(Loname, rs); st_matrixcolstripe(Loname, cs) st_matrix(Hiname, Hi); st_matrixrowstripe(Hiname, rs); st_matrixcolstripe(Hiname, cs) st_matrix(Qname, taus'); st_matrixcolstripe(Qname, cs) // posting matrices: e(b) is 1 x k row vector with var-name cols; // e(V) is k x k symmetric. We drop _cons from rnames // and append it at the end (Stata convention). string rowvector rn2 real matrix Vp, bp, perm real scalar i // _cons is rnames[1]; reorder to (vars..., _cons) rn2 = J(1, nb, "") for (i=2; i<=nb; i++) rn2[i-1] = rnames[i] rn2[nb] = "_cons" // permutation: position i in new layout came from position i+1 in old layout // for i 0) { Fsub = Fmat[(tstart-plag)..T, .] Fd = _nqar_build_F(Fsub, N, plag) } else Fd = J(rows(yv), 0, .) X = wylag, ylag, Zd, Fd kreg = cols(X) rnames = ("_cons", "WY", "WY_L1", "Y_L1") if (cols(Zmat) > 0) rnames = rnames, tokens(zlist) if (cols(Fmat) > 0) { fn = tokens(flist) for (k=0; k<=plag; k++) { for (j=1; j<=cols(Fmat); j++) { if (k==0) rnames = rnames, fn[j] else rnames = rnames, sprintf("%s_L%g", fn[j], k) } } } nb = cols(rnames) taus = strtoreal(tokens(qstr))' nq = rows(taus) user_grid = J(0, 1, .) gridmat = st_matrix(GRname) if (!(gridmat[1,1] == . & cols(gridmat)==1)) user_grid = vec(gridmat) B = J(nb, nq, .); SE = J(nb, nq, .) Tv = J(nb, nq, .); Pv = J(nb, nq, .) Lo = J(nb, nq, .); Hi = J(nb, nq, .) alphas = J(1, nq, .) gnorms = J(1, nq, .) zlev = invnormal(1 - (1-lev/100)/2) real matrix V_first real rowvector b_first V_first = J(nb, nb, .) b_first = J(1, nb, .) for (q=1; q<=nq; q++) { if (rows(user_grid) > 0) grid = user_grid else { Xpilot = wycur, X pilot = dnqr_qreg_via_st(yv, Xpilot, taus[q]) gam_init = pilot[1] pres = yv :- (J(rows(yv),1,1), Xpilot) * (pilot[cols(Xpilot)+1] \ pilot[1..cols(Xpilot)]') gam_se = max((0.01, sqrt(variance(pres))/sqrt(rows(yv)))) grid = gam_init :+ (range(-gpts, gpts, 1) :* (gscale*gam_se/gpts)) } found = dnqr_ivqr_grid(yv, wycur, X, Mz, grid, taus[q]) alpha = found[1] bfin = dnqr_qreg_via_st(yv :- alpha:*wycur, X, taus[q]) bcoef = J(nb, 1, .) bcoef[1] = bfin[kreg+1] bcoef[2] = alpha bcoef[3..nb] = bfin[1..kreg]' resid = yv :- alpha:*wycur :- (J(rows(yv),1,1), X) * (bcoef[1] \ bcoef[3..nb]) V = dnqr_powell_vcov(yv, wycur, X, Mz, resid, taus[q], btype, bscale) B[., q] = bcoef SE[., q] = sqrt(diagonal(V)) Tv[., q] = B[., q] :/ SE[., q] Pv[., q] = 2 :* (1 :- normal(abs(Tv[., q]))) Lo[., q] = B[., q] :- zlev :* SE[., q] Hi[., q] = B[., q] :+ zlev :* SE[., q] alphas[q] = alpha gnorms[q] = found[2] if (q==1) { b_first = bcoef' V_first = V } } netd = sum(W:>0) / (N*N - N) st_numscalar(Dname, netd) rs = J(nb, 1, ""), rnames' cs = J(nq, 1, ""), strofreal(taus, "%6.3f") st_matrix(Bname, B); st_matrixrowstripe(Bname, rs); st_matrixcolstripe(Bname, cs) st_matrix(SEname, SE); st_matrixrowstripe(SEname, rs); st_matrixcolstripe(SEname, cs) st_matrix(Tname, Tv); st_matrixrowstripe(Tname, rs); st_matrixcolstripe(Tname, cs) st_matrix(Pname, Pv); st_matrixrowstripe(Pname, rs); st_matrixcolstripe(Pname, cs) st_matrix(Loname, Lo); st_matrixrowstripe(Loname, rs); st_matrixcolstripe(Loname, cs) st_matrix(Hiname, Hi); st_matrixrowstripe(Hiname, rs); st_matrixcolstripe(Hiname, cs) st_matrix(Qname, taus'); st_matrixcolstripe(Qname, cs) st_matrix(Aname, alphas); st_matrixcolstripe(Aname, cs) st_matrix(Gname, gnorms); st_matrixcolstripe(Gname, cs) // post-able e(b)/e(V) -- reorder rnames so _cons is last string rowvector rn2 real matrix Vp, perm real rowvector bp real scalar i, j2 rn2 = J(1, nb, "") for (i=2; i<=nb; i++) rn2[i-1] = rnames[i] rn2[nb] = "_cons" perm = J(1, nb, .) for (i=1; i 1e-12) { S = luinv(I_n :- g1:*W) G = S * (g2:*W :+ g3:*I_n) } else G = g2:*W :+ g3:*I_n v0 = J(N, 1, 0) if (shocknode < 1 | shocknode > N) shocknode = 1 v0[shocknode] = shocksize IRF = dnqr_impulse(G, v0, H) nm = J(2, H+1, .) for (h=1; h<=H+1; h++) { nm[1, h] = sqrt(sum(IRF[., h]:^2)) nm[2, h] = max(abs(IRF[., h])) } st_matrix(IRFname, IRF) st_matrix(NORMname, nm) } // ===================================================================== // Simulator helpers // ===================================================================== real matrix _dnqr_W_dyad(real scalar N, real scalar delta) { real matrix A real scalar i, npairs, a, b, j A = J(N, N, 0) npairs = round(N*delta) for (i=1; i<=npairs; i++) { a = ceil(runiform(1,1)*N) b = ceil(runiform(1,1)*N) if (a==b) b = mod(b, N) + 1 A[a, b] = 1 A[b, a] = 1 } for (i=1; i<=N; i++) { if (sum(A[i, .])==0) { j = ceil(runiform(1,1)*N) if (j==i) j = mod(j, N) + 1 A[i, j] = 1 } } return(A) } real matrix _dnqr_W_powerlaw(real scalar N, real scalar alpha) { real matrix A real scalar i, nf, k real colvector pool, pick A = J(N, N, 0) for (i=1; i<=N; i++) { nf = ceil(runiform(1,1)^(-1/alpha)) nf = min((nf, N-1)) pool = (1::N) pool = select(pool, pool:!=i) _jumble(pool) pick = pool[1::min((nf, rows(pool)))] for (k=1; k<=rows(pick); k++) A[i, pick[k]] = 1 } for (i=1; i<=N; i++) { if (sum(A[i, .])==0) { k = mod(i, N) + 1 A[i, k] = 1 } } return(A) } real matrix _dnqr_W_block(real scalar N, real scalar K) { real matrix A real scalar i, j, blk, bi, bj, pr A = J(N, N, 0) blk = floor(N/K) if (blk < 1) blk = 1 for (i=1; i<=N; i++) { bi = floor((i-1)/blk) + 1 for (j=1; j<=N; j++) { if (i==j) continue bj = floor((j-1)/blk) + 1 pr = (bi==bj) ? 0.4 : 0.05 if (runiform(1,1) < pr) A[i, j] = 1 } } for (i=1; i<=N; i++) { if (sum(A[i, .])==0) { j = mod(i, N) + 1 A[i, j] = 1 } } return(A) } real matrix _dnqr_W_asym(real scalar N, real scalar h) { real matrix A real scalar i, j, hh A = J(N, N, 0) hh = round(h) for (i=1; i<=N; i++) { for (j=i+1; j<=min((N, i+hh)); j++) A[i, j] = 1 } for (i=1; i<=N; i++) { for (j=1; j 0) ? rnormal(N, nz, 0, 1) : J(N, 0, .) Fmat = (nf > 0) ? rnormal(T + burnin, nf, 0, 1) : J(T + burnin, 0, .) Tt = T + burnin Ymat = J(N, Tt, 0) for (t=2; t<=Tt; t++) { if (edist=="normal") err = rnormal(N, 1, 0, 1) else if (edist=="t") err = invnormal(uniform(N,1)) // approx; t is not built-in else err = rchi2(N, 1, edf) :- edf rhs = g2:*(W*Ymat[., t-1]) :+ g3:*Ymat[., t-1] :+ err if (nz > 0) rhs = rhs + Zmat * J(nz, 1, 0.3) if (nf > 0) rhs = rhs + J(N, 1, 1) :* (sum(Fmat[t, .]) * 0.2) if (abs(g1) > 1e-12) Ymat[., t] = luinv(I(N) :- g1:*W) * rhs else Ymat[., t] = rhs } Ymat = Ymat[., (burnin+1)..Tt] if (nf > 0) Fmat = Fmat[(burnin+1)..Tt, .] st_matrix(Wname, W) st_matrix(Yname, Ymat) if (nz > 0) st_matrix(Zname, Zmat) else st_matrix(Zname, J(0, 0, .)) if (nf > 0) st_matrix(Fname, Fmat) else st_matrix(Fname, J(0, 0, .)) } void _dnqr_sim_writeY(string scalar vn, real matrix Y) { real scalar N, T, i, t, row real colvector v N = rows(Y) T = cols(Y) v = J(N*T, 1, .) row = 1 for (t=1; t<=T; t++) { for (i=1; i<=N; i++) { v[row] = Y[i, t] row = row + 1 } } st_store(., vn, v) } void _dnqr_sim_writeZ(real matrix Z) { real scalar N, q, NT, k, i, t, row real colvector v if (rows(Z)==0 | cols(Z)==0) return N = rows(Z) q = cols(Z) NT = st_nobs() for (k=1; k<=q; k++) { v = J(NT, 1, .) row = 1 for (t=1; t<=NT/N; t++) { for (i=1; i<=N; i++) { v[row] = Z[i, k] row = row + 1 } } st_store(., sprintf("Z%g", k), v) } } void _dnqr_sim_writeF(real matrix F) { real scalar T, m, N, NT, k, i, t, row real colvector v if (rows(F)==0 | cols(F)==0) return T = rows(F) m = cols(F) NT = st_nobs() N = NT / T for (k=1; k<=m; k++) { v = J(NT, 1, .) row = 1 for (t=1; t<=T; t++) { for (i=1; i<=N; i++) { v[row] = F[t, k] row = row + 1 } } st_store(., sprintf("F%g", k), v) } } end