[Gretl-users] Fan chart
Someone asked me instructions on how to produce a "fan chart" (aka "rivers
of blood" chart) in gretl. These are typically used for density forecasts.
So I wrote a little function that I'm sharing here, if anyone is
interested. It's very rough, but I'm sure anyone with a decent knowledge
of hansl can customise it as per their needs.
Attention developers: I'0m not using the newly-introduced "band" option to
plot here, since you need more than one at the same time. Maybe we could
think about extending the syntax to "plot" in the future.
<hansl>
function string fanchart(series y, series se, matrix alphas)
matrix a = mreverse(sort(vec(alphas)))
scalar n = rows(a)
matrix X = {y , misszero(se)}
strings colcodes =
defarray("cccccc","888888","666666","333333")
scalar rand = ceil(1000*muniform(1,1))
string gfname = sprintf("\"(a)dotdir/tmp%03d.gp\"", rand)
string dfname = sprintf("(a)dotdir/tmp%03d.dat", rand)
outfile @gfname --write --quiet
printf "set nokey\nset ylabel '%s'\n", argname(y)
printf "plot \\\n"
loop i = 1..n --quiet
z = invcdf(z, 1 - (1-a[i])/2)
string code = colcodes[i]
printf "'%s' using 1:($2-%g*$3):($2+%g*$3) notitle lc rgb
\"0x%s\" w filledcurve, \\\n", \
dfname, z, z, code
endloop
printf "'%s' using 1:2 w lines lt 1\n", dfname
outfile --close
outfile @dfname --write --quiet
loop i = 1..rows(X) --quiet
scalar tt = $obsmajor[i] + ($pd>1 ? ($obsminor[i] - 1)/$pd : 0)
printf "%g\t%g\t%g\n", tt, X[i,1], X[i,2]
endloop
outfile --close
return gfname
end function
# usage example
nulldata 48
setobs 12 1990:1 --special-time-series
scalar n = $nobs
series y = 0.1*time + filter(0.5*mnormal(n,1), 1, {1.8, -0.9})
series se = cum(time<$nobs/1.5 ? NA : uniform())
matrix alphas = {.83, 0.95, 0.99}
string haha = fanchart(y, se, alphas)
gnuplot --input=@haha --output=display
</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
-------------------------------------------------------