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>

Thanky so much for your answers.
Currently I'm doing J-N plots on excel, but the time spent and the probability of errors is very high, so it would be interesting having Gretl doing the hard work.
This author http://mattgolder.com/interactions provides Stata code for ploting J-N, both for a simple interaction and for three-way interactions as well.

Best,
Marcelo