[Gretl-users] Re: Johnson-Neyman plots
On Tue, 7 Mar 2023, Riccardo (Jack) Lucchetti wrote:
One of the things that would make this easier to do in hansl would be
the
possibility (now absent, but discussed some time ago in a developer meeting
IIRC) to have more than one "band" specification in the gnuplot/plot
command.
Anyway: if anyone's interested in developing a package for this kind of
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 = sd(x)
x0 = min(x) - 0.025*sdx
x1 = max(x) + 0.025*sdx
matrix X = x0 + seq(0, n-1) * ((x1 -x0)/(n-1))
X = 1 ~ X'
matrix m = X * coef
matrix v = 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 = mod.coeff[select]
matrix V = mod.vcv[select, select]
matrix M = meff_mat(coef, V, x)
PLT = M[,{2,1}]
B = M[,{2,3}]
plot PLT
option with-lines
options band=B,1.96 band-style=fill
end plot --output=display
end function
### example
set verbose off
set seed 11080
nulldata 96
x = normal()
z = normal()
xz = x * z
y = 1 + x - z + 0.25*xz + 0.8 * normal()
ols y const x z xz
sel = {2,4}
JNplot(z, $model, sel)
</hansl>
-------------------------------------------------------
Riccardo (Jack) Lucchetti
Dipartimento di Scienze Economiche e Sociali (DiSES)
Università Politecnica delle Marche
(formerly known as Università di Ancona)
r.lucchetti(a)univpm.it
http://www2.econ.univpm.it/servizi/hpp/lucchetti
-------------------------------------------------------