$title Use GDXRANK to report multisectoral Monte Carlo results (rank03,SEQ=137)
$onText
One of the most perplexing challenges in economic modeling with GAMS is to
present multisectoral results in an easily interpreted format. One simple idea
is to present sectoral results in a sorted sequence to make it easier to
identify the most seriously affected sectors. The presentation of results in a
multisectoral model is made even more challenging when model results are
generated for a randomized set of scenarios. A summary of Monte Carlo results
involves reporting both mean results and their sensitivity. One means of
characterizing the sensitivity of model results is to report functions of the
sample distribution such as the upper and lower quartiles.
This example illustrates how $libInclude rank can be used to help report results from
the Monte Carlo analysis of a multisectoral model.
Keywords: rank, Monte-Carlo, quartiles
$offText
Set
run 'Samples (potentially not all solved)' / 1*1000 /
i 'Sectors for which results are to be compared'
/ELE,OLE,OLP,GAS,COA,OFU,FME,NFM,CHM,MWO,TPP
CNM,CLO,FOO,OTH,CON,AGF,RLW,TRK,PIP,MAR,AIR
TRO,TMS,PST,TRD,CAT,OIN,PSM,SSM,ECM,SCS,GEO,FIN,ADM /;
Parameter v(run,i) 'Monte-Carlo results for all sectors';
* Load the sectoral results from the sensitivity analysis from
* a GDX file:
$gdxIn 'ex3.gdx'
$load v
Set
s(run) 'Solved cases'
k 'Count used for quartiles' / 1*1000 /
ki 'Count used for sectors' / 1*100 /
qtl 'Quartiles (set)' / q25, q50, q75 /
imap(ki,i) 'Mapping from ordered list to sector labels';
Parameter
qvalue(i,*) 'Quartile values'
mean(i) 'Mean impact on sector i'
meanrank(i) 'Mean rank of sector i'
x(run) 'Vector used for sorting'
r(run) 'Rank values returned'
quartile(qtl) 'Quartiles (evaluated)'
qv(qtl) 'Quartiles (values)' / q25 25, q50 50, q75 75 /;
* Identify which scenarios have been solved (some of the runs
* may have failed):
s(run) = yes$(smax(i,abs(v(run,i))) <> 0);
* Evaluate means and support for results:
qvalue(i,"mean") = sum(s, v(s,i))/card(s);
qvalue(i,"min") = smin(s, v(s,i));
qvalue(i,"max") = smax(s, v(s,i));
display qvalue;
* Determine ranking of sectors by mean impact:
mean(i) = qvalue(i,"mean");
$libInclude rank mean i meanrank
* The following statement creates a tuple matching the ordered
* set, ki, to the set of sectors, i. In this tuple, the sequence of
* assignments corresponds to increasing mean impacts:
imap(ki+(meanrank(i)-ord(ki)),i) = yes;
* Evaluate quartiles of sectoral impacts for each sector:
loop(i,
x(s) = v(s,i);
* Load quartile with the perctiles to be
* evaluated (25th, 50th and 75th):
quartile(qtl) = qv(qtl);
$ libInclude rank x s r quartile
* Save the quartile values:
qvalue(i,qtl) = quartile(qtl);
);
display qvalue;
Parameter results(ki,i,*) 'Summary of impacts (sorted)';
results(ki,i,"mean")$imap(ki,i) = mean(i);
results(ki,i,qtl)$imap(ki,i) = qvalue(i,qtl);
display results;
* Use existence of plot.gms in inclib in the GAMS system directory as indicator for installed gnuplot
$if not exist "%gams.sysdir%inclib/plot.gms" $exit
* Write out a GNUPLOT file to generate a chart of the results:
File kplt / ex3.gnu /;
put kplt;
kplt.lw = 0;
put "reset"/;
put 'set title "Sectoral Impacts with Quartiles"'/;
put "set linestyle 1 lt 8 lw 1 pt 8 ps 0.5"/;
put "set grid"/;
put 'set ylabel "% change"'/;
put "set xzeroaxis"/;
put "set bmargin 4"/;
put "set xlabel 'sector'"/;
put 'set xrange [1:',card(i),']'/;
put 'set xtics rotate (';
loop(ki, loop(i$imap(ki,i), put '\'/' "',i.tl,'" ', ord(ki):0,',';));
put @(file.cc-1) ')'/;
put "plot 'ex3.dat' notitle with errorbars ls 1"/;
putClose;
File kpltdata / ex3.dat /;
put kpltdata;
kpltdata.nr = 2;
kpltdata.nw = 14;
kpltdata.nd = 6;
loop(ki, loop(i$imap(ki,i), put ord(ki):0,qvalue(i,"mean"), qvalue(i,"q25"), qvalue(i,"q75")/;));
putClose;
execute 'wgnupl32 ex3.gnu -';