*! version 1.0.0 11jul2026

Title

hpcm methods -- the algorithm behind hpcm, equation by equation

Overview

Let w(t)=(x(t),y(t),z(t)) be a zero-mean second-order stationary process,
p1+p2+p3 = K dimensional. hpcm fits a VAR(p)

A(L) w(t) = e(t), Cov(e(t)) = Sy, A(L)=I-A1 L-...-Ap L^p

by OLS on the centered data, and works entirely from A(L) and Sy. Index the
blocks o=(x,y) [size q=p1+p2] and 3=z [size p3]. All spectra use the convention
f({&lambda}) = (1/2{&pi}) {&Psi}(e^-i{&lambda}) Sy {&Psi}(e^-i{&lambda})*
with {&Psi}(z)=A(z)^-1, so the innovation-representation transfer function is
{&Lambda}(z)=A(z)^-1 Sy^1/2 and {&Lambda}(0){&Lambda}(0)'=Sy.

The one-way-effect component of z

z0,0,-1(t) is the residual of z(t) projected on H.
At the innovation level this is simply the z-innovation orthogonalized on the
contemporaneous (x,y)-innovations:

z0,0,-1(t) = P e(t), P = [ -Sy_3o Sy_oo^-1 , I(p3) ]

a K-selecting p3xK matrix. It is white noise with covariance
Sy_33.o = Sy_33 - Sy_3o Sy_oo^-1 Sy_o3 (the Schur complement), matching the
paper's ~f33({&lambda}) = (1/2{&pi}) Sy_33.o.

Cleaning x and y (eq 4.2)

u(t) and v(t) are the residuals of x(t) and y(t) after projecting out all
leads and lags of z0,0,-1. Removing the P-direction from the innovation covariance
gives the annihilated covariance

M Sy = Sy - Sy P' (Sy_33.o)^-1 P Sy (rank q),

and the cleaned pair is (u,v)(t) = {&Psi}_o(L) M e(t), whose spectral density is

g({&lambda}) = (1/2{&pi}) {&Psi}_o(e^-i{&lambda}) (M Sy) {&Psi}_o(e^-i{&lambda})*,

with {&Psi}_o = the first q rows of A(e^-i{&lambda})^-1. This is algebraically
identical to Hosoya's eq (4.2), g = f_oo - 2{&pi} ~f_o3 ~f33^-1 ~f_3o, but is
computed in the numerically clean annihilated form above (the two agree because
M Sy M' = M Sy).

The partial measures

Partition g into u-block (1) and v-block (2). With {&Sigma} the one-step
prediction-error covariance of (u,v) and B(e^-i{&lambda})={&Gamma}(0)
{&Gamma}(e^-i{&lambda})^-1 its canonical factor (below):

Directional (eq 4.5), with Q_b=[-{&Sigma}21{&Sigma}11^-1, I],
{&Sigma}22:1={&Sigma}22-{&Sigma}21{&Sigma}11^-1{&Sigma}12, and
~g21({&lambda})=Q_b B(e^-i{&lambda}) g({&lambda})[:,1:p1]:

PM(y{&rarr}x:z)({&lambda}) = ln det g11 - ln det[ g11 - 2{&pi} ~g12 {&Sigma}22:1^-1 ~g21 ]

PM(x{&rarr}y:z)({&lambda}) is the mirror image (swap blocks 1 and 2).

Association (eq 4.6):

PM(x,y:z)({&lambda}) = ln det g11 + ln det g22 - ln det g

Reciprocity (Theorem 4.2, by identity):

PM(x:y:z)({&lambda}) = PM(x,y:z)({&lambda}) - PM(x{&rarr}y:z)({&lambda}) - PM(y{&rarr}x:z)({&lambda})

Each overall measure is (1/{&pi}){&int}_0^{&pi} M({&lambda})d{&lambda} (the
integrand is even in {&lambda}; Theorem 3.1 / 4.1), evaluated by the trapezoidal
rule on the grid, augmented with the band() endpoints so band integrals are
exact at the edges. All measures are floored at 0 (they are non-negative in theory;
tiny negative values from quadrature are set to 0).

Canonical factorization of the cleaned pair

The directional measures need {&Gamma}(0){&Gamma}(e^-i{&lambda})^-1. For
a finite VAR this equals the lag polynomial evaluated on the unit circle, so
hpcm represents (u,v) by a long VAR(m): it forms the autocovariances

R(h) = {&Sigma}_j {&Psi}_o(j+h) (M Sy) {&Psi}_o(j)', {&Psi}(j)=companion MA weights,

solves the multivariate Yule-Walker system Brow = rvec (block-Toeplitz)^-1 for
B1..Bm and {&Sigma}=R(0)-Brow rvec', and sets B(e^-i{&lambda})=I-{&Sigma}Bk
e^-ik{&lambda}. As m grows this converges (geometrically for a stable VAR) to
the exact canonical factor. mlag() sets m; matrunc() truncates the R(h)
sum.

Wald inference (Section 6)

Let {&theta}=(vec A1..Ap, vech Sy) with finite-sample covariance V (OLS:
V_A={&Sigma}{&otimes}(Z'Z)^-1 for the lag block, V_S=(1/T) form for vech Sy,
independent under Gaussianity). Let G({&theta}) collect the band measures
[PM(y{&rarr}x:z), PM(x{&rarr}y:z)]. With the numerical Jacobian
D{&theta}G (central differences) the paper's statistic reduces to

W = G' [ D{&theta}G V D{&theta}G' ]^-1 G ~ {&chi}²(m),

(the factor T of the textbook form cancels the T in V). hpcm reports the two
single-direction tests (m=1) and the joint test (m=2). Under H0 a measure lies on
the boundary 0, so D{&theta}G may be near zero and the {&chi}² is conservative;
this is a documented property of the method.

Step -> equation map

code step paper reference

hpcm_fit (OLS VAR of x,y,z) Section 5, eq (5.5); Sec 3 factorization

hpcm_msig (P, Sy_33.o, M Sy) Sec 2-4, one-way component; ~f33

g = (1/2{&pi}) {&Psi}_o M Sy {&Psi}_o* eq (4.2)

hpcm_acov + hpcm_yw Sec 3 canonical factorization; Hosoya (1991)

directional M({&lambda}) eq (3.9)-(3.11), (4.4)-(4.5)

association M({&lambda}) eq (4.6)

reciprocity M({&lambda}) Theorem 4.2, eq (4.8)

overall = (1/{&pi}){&int}M Theorem 3.1 / 4.1

difference() Sec 5 reproducible reduction, eq (5.4)

Wald test Section 6

Author

Dr Merwan Roudane

merwanroudane920@gmail.com

github.com/merwanroudane