open denmark.gdt
list LST = LRM LRY IBO IDE
scalar p = 2 # lags in VAR
scalar USE_F = 1 # small sample correction?

var p LST --quiet

k = nelem(LST)
matrix theta = vec($coeff)
matrix V = $sigma ** $xtxinv
if USE_F
  scalar df = $T - $ncoeff
  V *= $T/df
endif

matrix GC = zeros(k, k)
cnameset(GC, LST)
rnameset(GC, LST)

matrix idx = seq(1,p) + 1
loop i = 1..k
  loop j = 1..k
    GC[i,j] = qform(theta[idx]', invpd(V[idx,idx]))
    idx += (j==k)? p+1 : p
  endloop
endloop

if USE_F
  GC /= p
  matrix pvals = pvalue(F, p, df, GC)
else
  matrix pvals = pvalue(X, p, GC)
endif

cnameset(pvals, LST)
rnameset(pvals, LST)
print GC pvals