On Tue, 7 Mar 2023, Riccardo (Jack) Lucchetti wrote:

> One of the things that would make this easier to do in hansl would be t=
he=20
> possibility (now absent, but discussed some time ago in a developer mee=
ting=20
> IIRC)  to have more than one "band" specification in the gnuplot/plot=20
> command.

Anyway: if anyone's interested in developing a package for this kind of=20
plots, here's a rough attempt to a starting point.

<hansl>
function matrix meff_mat(matrix coef, matrix vcv,
                          const series x, scalar n[32])

     scalar sdx =3D sd(x)
     x0 =3D min(x) - 0.025*sdx
     x1 =3D max(x) + 0.025*sdx
     matrix X =3D x0 + seq(0, n-1) * ((x1 -x0)/(n-1))
     X =3D 1 ~ X'
     matrix m =3D X * coef
     matrix v =3D sumr(X .* (X*vcv))
     return X[,2] ~ m ~ sqrt(v)
end function

function void JNplot(const series x, const bundle mod, matrix select,
                      scalar n[32])

     matrix coef =3D mod.coeff[select]
     matrix V =3D mod.vcv[select, select]
     matrix M =3D meff_mat(coef, V, x)

     PLT =3D M[,{2,1}]
     B =3D M[,{2,3}]

     plot PLT
         option with-lines
         options band=3DB,1.96 band-style=3Dfill
     end plot --output=3Ddisplay

end function

### example

set verbose off
set seed 11080
nulldata 96

x =3D normal()
z =3D normal()
xz =3D x * z
y =3D 1 + x - z + 0.25*xz + 0.8 * normal()

ols y const x z xz
sel =3D {2,4}
JNplot(z, $model, sel)
</hansl>