$ontext
    analyseSample.gms
    File for analysing results of estimations carried out with program
    generateSampleGrid.gms and stored in GDX-files smp1, smp2, ... smpN

    IMPORTANT: set the scalar "npr" to the number of processes that were
    run in parallel with generateSampleGrid!

    Torbjorn Jansson
    LEI, The Hague, NL
$offtext
$offlisting
$setglobal PAPEROUT ON
$setglobal POSTEROUT OFF
$set TRAD a8
$set BLEP a10
$set OBS  a0
$set repdir rep
scalar npr "Number of gdx-files to load results from" /4/;

file con /con/; con.lw = 0; con.tw = 0; con.nw = 0; con.nd = 0;
file batch /titlebatch.bat/; batch.lw = 0;

set m  Models generated ; alias(m,m1);
set n  Estimation attempts per model;
set a  Solution approaches available in total; alias(a,a1,a2,a3);
set i  Regions ;
set t  Periods;

alias(i,j,k,l);
alias(t,t1);

putclose con / "... Reading sets ..." /;

$gdxin "smp1.gdx"
$load m n a i t
$gdxin

set aon(a) Solution approaches to analyse /a8,a10/;

parameter ptru(m,t,i)     True price;
parameter ctru(m,i,j)     True transportation cost;
parameter stru(m,i)       True storage cost;
parameter etru(m,t,i)     True excess demand;

parameter pest(m,n,t,i,a) Estimated price;
parameter cest(m,n,i,j,a) Estimated transportation cost;
parameter sest(m,n,i,a)   Estimated storage cost;

variable VDUM; equation EDUM; EDUM .. VDUM =e= 10;
model MDUM 'Dummy model clears memory for some odd reason' /EDUM/;
MDUM.solprint = 2;

set mrun(m) Particular model to analyse; mrun(m) = no;

putclose con / "... Reading data ..." /;

execute_load "smp1.gdx" ptru ctru stru etru pest cest sest;

parameter pest1(m,n,t,i,a) Estimated price temporary;
parameter cest1(m,n,i,j,a) Estimated transportation cost temporary;
parameter sest1(m,n,i,a)   Estimated storage cost temporary;

* Load and merge data from up to four different gdx-files
if(npr>1,
    execute_load "smp2.gdx" mrun pest1=pest cest1=cest sest1=sest;
    pest(mrun,n,t,i,a) = pest1(mrun,n,t,i,a);
    cest(mrun,n,i,j,a) = cest1(mrun,n,i,j,a);
    sest(mrun,n,i,a)   = sest1(mrun,n,i,a);
    option kill=pest1; option kill=cest1; option kill=sest1; mrun(m) = no;
);

if(npr>2,
    execute_load "smp3.gdx" mrun pest1=pest cest1=cest sest1=sest;
    pest(mrun,n,t,i,a) = pest1(mrun,n,t,i,a);
    cest(mrun,n,i,j,a) = cest1(mrun,n,i,j,a);
    sest(mrun,n,i,a)   = sest1(mrun,n,i,a);
    option kill=pest1; option kill=cest1; option kill=sest1; mrun(m) = no;
);

if(npr>3,
    execute_load "smp4.gdx" mrun pest1=pest cest1=cest sest1=sest;
    pest(mrun,n,t,i,a) = pest1(mrun,n,t,i,a);
    cest(mrun,n,i,j,a) = cest1(mrun,n,i,j,a);
    sest(mrun,n,i,a)   = sest1(mrun,n,i,a);
    option kill=pest1; option kill=cest1; option kill=sest1; mrun(m) = no;
);

*execute_unload "smpt.gdx" pest cest sest m i t a;
*$exit

scalar count /0/;
scalar np Number of prices estimated;
scalar nc Number of transport costs estimated;
scalar ns Number of storage costs estimated;
scalar nn Number of iterations per model;

np = card(i)*card(t);
nc = card(i)*(card(i)-1)/2;
ns = card(i);
nn = card(n);

set AD(i,j) Admissible trade flows;
AD(i,j) = yes $ (ord(j) gt ord(i));

parameter cord(i,j) Ordering of flows;

count=0;
loop(AD, count=count+1; cord(AD)=count);

scalar son Storage on /0/;
son = card(t)-1;

set obs/o1*o10000/;
set uobs(obs);

* Some basic reporting to list file

option count:0;
count = card(m);
display "Number of true models in sample:",count;
count = card(n);
display "Number of solutions per model:",count;
count = card(t);
display "Number of time periods in storage model (1=no storage)",count;
count = card(i);
display "Number of regions in model",count;
count = sum((m,n,i,j) $ (AD(i,j) and (1 = cest(m,n,i,j,"%OBS%"))),1);
display "Number of truncated cost items:", count;
count = sum((m,n,t,i) $ ((139 = pest(m,n,t,i,"%OBS%"))
                      or (101 = pest(m,n,t,i,"%OBS%"))),1);
display "Number of truncated price items:", count;

* Select methods to analyse and discard data for other methods

pest(m,n,t,i,a) $ (not aon(a)) = 0;
cest(m,n,i,j,a) $ (not aon(a)) = 0;
sest(m,n,i,a)   $ (not aon(a)) = 0;
solve MDUM using NLP maximising VDUM;
set compare(a,a1) Pairs of methods to compare;
compare(a,a1) $ (aon(a) and aon(a1) and (not sameas(a,a1))) = yes;
set compnow(a,a1) Pairt of methods currently compared;
    compnow(a,a1) = no;


* ----------------------------------------------------------------------
*   Declare statistical measures to compute
* ----------------------------------------------------------------------

set parvec /preg Regional price, tpc Transport cost, stc Storage cost/;
set meth /TRAD Traditional method , BLEP Bilevel Estimation Program/;

parameter pmean(m,t,i,a)    Sample mean;
parameter cmean(m,i,j,a);
parameter smean(m,i,a);

parameter psvar(m,t,i,a)    Sample variance;
parameter csvar(m,i,j,a);
parameter ssvar(m,i,a);

parameter ptval(m,t,i,a)    T-values for test estimated equals true;
parameter ctval(m,i,j,a);
parameter stval(m,i,a);

set alpha 'Levels of significanse' /alpha950,alpha990,alpha999/;
parameter pcrit(alpha) 'Critical test statistic value' /
    alpha950    0.05
    alpha990    0.01
    alpha999    0.001 /;

parameter tsummary(a,*);

set mse Items used for analysis of mean squared error /
    msetot  Mean squared error
    msebias Bias
    msevar  Variance    /;

set totmse Items used for analysis of mean squared error /
    MMSE    Mean MSE
    MSBIAS  Mean Squared Bias
    MVAR    Mean (pooled) Variance    /;

parameter pmse(m,t,i,a,mse)  Mean squared error per item;
parameter cmse(m,i,j,a,mse);
parameter smse(m,i,a,mse);

parameter psrmse(*,a,totmse)    Sum of root mean squared error;
parameter csrmse(*,a,totmse);
parameter ssrmse(*,a,totmse);

parameter pmmse(m,t,i,t,i,a)    "Full MSE matrix, price block";
parameter cmmse(m,i,i,i,i,a)    "Full MSE matrix, cost block";
parameter pcmmse(m,t,i,i,i,a)   "Full MSE matrix, price-cost block";

option pmse:3:3:2;
option cmse:3:3:2;
option smse:3:2:2;

option psrmse:3:1:2;
option csrmse:3:1:2;
option ssrmse:3:1:2;

parameter bartlett(m,a,parvec) Bartlett statistic equality of variances;
parameter levene(m,a,parvec)   Levene statistic equality of variances;

option bartlett:2:1:2;
option levene:2:1:2;

*   Items for GDX-rank

parameter pct(*) /median 50.0/;
parameter someSort1(*);
parameter someSort2(*);
parameter someRank(*);

* ----------------------------------------------------------------------
*   Compute mean and sample variance of estimated parameters
* ----------------------------------------------------------------------

putclose con / "... computing means ..." /;

pmean(m,t,i,aon)           = sum(n, pest(m,n,t,i,aon))/card(n);
cmean(m,i,j,aon) $ AD(i,j) = sum(n, cest(m,n,i,j,aon))/card(n);
smean(m,i,aon)   $ son     = sum(n, sest(m,n,i,aon))  /card(n);

putclose con / "... computing variances ..." /;

psvar(m,t,i,aon)           = sum(n, sqr(pest(m,n,t,i,aon)
                                        - pmean(m,t,i,aon)))/(card(n)-1);
csvar(m,i,j,aon) $ AD(i,j) = sum(n, sqr(cest(m,n,i,j,aon)
                                        - cmean(m,i,j,aon)))/(card(n)-1);
ssvar(m,i,aon)   $ son     = sum(n, sqr(sest(m,n,i,aon)
                                        - smean(m,i,aon)))  /(card(n)-1);

putclose con / "... computing t-values ..." /;

ptval(m,t,i,aon)            = (ptru(m,t,i)-pmean(m,t,i,aon))
                              /sqrt(psvar(m,t,i,aon)/card(n));
ctval(m,i,j,aon) $ AD(i,j)  = (ctru(m,i,j)-cmean(m,i,j,aon))
                              /sqrt(csvar(m,i,j,aon)/card(n));
stval(m,i,aon)   $ son      = (stru(m,i)  -smean(m,i,aon))
                              /sqrt(ssvar(m,i,aon)/card(n));

*execute_unload "test.gdx";
*$exit

* ----------------------------------------------------------------------
*   Compute mean squared error and bias of estimators
* ----------------------------------------------------------------------

putclose con / "... computing MSE ..." /;

pmse(m,t,i,aon,"msebias")           = pmean(m,t,i,aon)-ptru(m,t,i);
cmse(m,i,j,aon,"msebias") $ AD(i,j) = cmean(m,i,j,aon)-ctru(m,i,j);
smse(m,i,aon,"msebias")   $ son     = smean(m,i,aon)  -stru(m,i);

*option kill = ptru; option kill = ctru;

pmse(m,t,i,aon,"msevar")            = psvar(m,t,i,aon);
cmse(m,i,j,aon,"msevar") $ AD(i,j)  = csvar(m,i,j,aon);
smse(m,i,aon,"msevar")   $ son      = ssvar(m,i,aon);

option kill = psvar; option kill = csvar; option kill = ssvar;

pmse(m,t,i,aon,"msetot")            = pmse(m,t,i,aon,"msevar")
                                    + sqr(pmse(m,t,i,aon,"msebias")) ;
cmse(m,i,j,aon,"msetot") $ AD(i,j)  = cmse(m,i,j,aon,"msevar")
                                    + sqr(cmse(m,i,j,aon,"msebias")) ;
smse(m,i,aon,"msetot")   $ son      = smse(m,i,aon,"msevar")
                                    + sqr(smse(m,i,aon,"msebias"))   ;

putclose con / "... computing sum of root(MSE) ..." /;

psrmse(m,aon,"msbias") = sum((t,i), sqr(pmse(m,t,i,aon,"msebias")))
                        /(card(i)*card(t));
csrmse(m,aon,"msbias") = sum((i,j)$ AD(i,j),sqr(cmse(m,i,j,aon,"msebias")))
                        /(card(i)*(card(i)-1)/2);
ssrmse(m,aon,"msbias") = sum(i$ son, sqr(smse(m,i,aon,"msebias")))
                        /card(i);;

*   Compute pooled variance

psrmse(m,aon,"mvar") = sum((t,i),           pmse(m,t,i,aon,"msevar"))
                       /(card(i)* card(t));
csrmse(m,aon,"mvar") = sum((i,j) $ AD(i,j), cmse(m,i,j,aon,"msevar"))
                       /(card(i)*(card(i)-1)/2);
ssrmse(m,aon,"mvar") = sum( i    $ son,     smse(m,i,aon,"msevar"))
                       /card(i);

*   Compute sum of Root Mean Squared Error, SRMSE

psrmse(m,aon,"mmse")  = sum((t,i),           pmse(m,t,i,aon,"msetot"))
                        /(card(i)* card(t));
csrmse(m,aon,"mmse")  = sum((i,j) $ AD(i,j), cmse(m,i,j,aon,"msetot"))
                        /(card(i)*(card(i)-1)/2);
ssrmse(m,aon,"mmse")  = sum( i    $ son,     smse(m,i,aon,"msetot"))
                        /card(i);

*   Average over all models

psrmse("tot",aon,totmse)   = sum(m, psrmse(m,aon,totmse));
csrmse("tot",aon,totmse)   = sum(m, csrmse(m,aon,totmse));
ssrmse("tot",aon,totmse)   = sum(m, ssrmse(m,aon,totmse));

putclose con / "Start analysing MSE superiority:" /;

*   Compute full MSE matrix

loop(aon,
    putclose con / "... of prices, ",aon.tl /;
    pmmse(m,t,i,t1,j,aon) $ (ord(i) le ord(j))
        = sum(n, (pest(m,n,t ,i,aon) - ptru(m,t ,i))
                *(pest(m,n,t1,j,aon) - ptru(m,t1,j)))
        / card(n);

    putclose con / "... of costs, ",aon.tl /;
    cmmse(m,i,j,k,l,aon) $ ((AD(i,j) and AD(k,l)) and (cord(i,j) le cord(k,l)))
        = sum(n, (cest(m,n,i,j,aon) - ctru(m,i,j))
                *(cest(m,n,k,l,aon) - ctru(m,k,l)))
        / card(n);

    putclose con / "... of price-costs interaction, ",aon.tl /;
    pcmmse(m,t,i,j,k,aon) $ AD(i,j)
        = sum(n, (pest(m,n,t ,i,aon) - ptru(m,t ,i))
                *(cest(m,n,j,k,aon)  - ctru(m,j,k)))
        / card(n);
);

*   Check for strong MSE superiority by checking if MSE(a1)-MSE(a2) is nsd.

set supset /strong, weak/;
set partit /p,c,pc/;
parameter msesup(*,a,a1,supset,partit)
variable xp(t,i), xc(i,j), z;
equation eqf Quadratic form using partitioned MSE matrix;
equation erz Constraint on objective value;

eqf.. z =E=
    sum((m,a,a1) $ (mrun(m) and compnow(a,a1)),
      sum((t,i,t1,j) $ (ord(i) le ord(j)),
        xp(t,i)*(pmmse(m,t,i,t1,j,a)-pmmse(m,t,i,t1,j,a1))*xp(t1,j))
   +  sum((t,i,t1,j) $ (ord(i) gt ord(j)),
        xp(t,j)*(pmmse(m,t,j,t1,i,a)-pmmse(m,t,j,t1,i,a1))*xp(t1,i))
   +  sum((i,j,k,l) $ ((AD(i,j) and AD(k,l)) and (cord(i,j) le cord(k,l))),
        xc(i,j)*(cmmse(m,i,j,k,l,a)-cmmse(m,i,j,k,l,a1))*xc(k,l))
   +  sum((i,j,k,l) $ ((AD(i,j) and AD(k,l)) and (cord(i,j) gt cord(k,l))),
        xc(k,l)*(cmmse(m,k,l,i,j,a)-cmmse(m,k,l,i,j,a1))*xc(i,j))
   +2*sum((t,i,j,k),
        xp(t,i)*(pcmmse(m,t,i,j,k,a)-pcmmse(m,t,i,j,k,a1))*xc(j,k))
    );

erz.. z =L= 1000;

model checknsd /eqf,erz/;
    checknsd.solprint  = 2;
    checknsd.solvelink = 2;
    checknsd.optfile = 2;

$onecho > conopt.op2
lfilos =  100
lfilog =  100
$offecho

loop(compare(a2,a3),
    loop(m1,
        putclose batch "title Superiority: ",m1.tl,", check ",a2.tl,"-",a3.tl;
        execute "titlebatch";

        mrun(m1) = yes;
        compnow(a2,a3) = yes;

*       Analyse price block for strong MSE superiority
        xp.lo(t,i)=-inf; xp.up(t,i)=inf;
        xp.l(t,i)=1; xc.fx(i,j)=0;
        solve checknsd using nlp maximising z;
        msesup(m1,a2,a3,"strong","p") $ (z.l le 500) = yes;

*       Analyse cost block for strong MSE superiority
        xc.lo(i,j)=-inf;xc.up(i,j)=inf;
        xp.fx(t,i)=0; xc.l(i,j)=1;
        solve checknsd using nlp maximising z;
        msesup(m1,a2,a3,"strong","c") $ (z.l le 500) = yes;

*       Analyse full matrix for strong MSE superiority
        if(msesup(m1,a2,a3,"strong","p") and msesup(m1,a2,a3,"strong","c"),
            xc.lo(i,j)=-inf;xc.up(i,j)=inf;
            xp.lo(t,i)=-inf; xp.up(t,i)=inf;
            xp.l(t,i)=1;xc.l(i,j)=1;
            solve checknsd using nlp maximising z;
            msesup(m1,a2,a3,"strong","pc") $ (z.l le 500) = yes;
        );

*       Analyse price block for weak MSE superiority
        z.l = sum((t,i), pmmse(m1,t,i,t,i,a2)-pmmse(m1,t,i,t,i,a3));
        msesup(m1,a2,a3,"weak","p") $ (z.l le 0) = yes;

*       Analyse cost block for weak MSE superiority
        count = z.l;
        z.l = sum((i,j) $ AD(i,j), cmmse(m1,i,j,i,j,a2)-cmmse(m1,i,j,i,j,a3));
        msesup(m1,a2,a3,"weak","c") $ (z.l le 0) = yes;

*       Analyse full matrix for weak MSE superiority
        z.l = count + z.l;
        msesup(m1,a2,a3,"weak","pc") $ (z.l le 0) = yes;

*       if(msesup(m1,a2,a3,"weak","p") and msesup(m1,a2,a3,"weak","c"),
*           msesup(m1,a2,a3,"weak","pc") $ (z.l le 0) = yes;
*       );

        compnow(a2,a3) = no;
        mrun(m1) = no;
    );
);

msesup("tot",a,a1,supset,partit) = sum(m, msesup(m,a,a1,supset,partit));
execute_unload "mse.gdx" pmmse,cmmse,pcmmse,msesup;
msesup(m,a,a1,supset,partit) = 0;

display "Pair-wise analyses of MSE superiority:", msesup;

*   Clean up memory

option kill = ptru;  option kill = ctru;
option kill = pmmse; option kill = cmmse;
option kill = pcmmse;option kill = msesup;

*$exit
* ----------------------------------------------------------------------
*   Declare some parameters to hold values to plot
* ----------------------------------------------------------------------

parameter plotpmse(m,meth)  Plot of SRMSE for prices for Poster;
parameter plotcmse(m,meth)  Plot of SRMSE for costs for Poster;

parameter plotpvar(m,meth)    Plot of pooled price variance for Poster;
parameter plotcvar(m,meth)    Plot of pooled price variance for Poster;

parameter plotpbias(m,meth)    Plot of price ABIAS for Poster;
parameter plotcbias(m,meth)    Plot of cost ABIAS for Poster;

parameter plotphomo(obs,*)  Plot of individual price variance estimates;
parameter plotchomo(obs,*)  Plot of individual cost variance estimates;

plotpmse(m,"TRAD") = psrmse(m,"%TRAD%","mmse");
plotcmse(m,"TRAD") = csrmse(m,"%TRAD%","mmse");
plotpmse(m,"BLEP") = psrmse(m,"%BLEP%","mmse");
plotcmse(m,"BLEP") = csrmse(m,"%BLEP%","mmse");

plotpvar(m,"TRAD") = psrmse(m,"%TRAD%","mvar");
plotcvar(m,"TRAD") = csrmse(m,"%TRAD%","mvar");
plotpvar(m,"BLEP") = psrmse(m,"%BLEP%","mvar");
plotcvar(m,"BLEP") = csrmse(m,"%BLEP%","mvar");

plotpbias(m,"TRAD") = psrmse(m,"%TRAD%","msbias");
plotcbias(m,"TRAD") = csrmse(m,"%TRAD%","msbias");
plotpbias(m,"BLEP") = psrmse(m,"%BLEP%","msbias");
plotcbias(m,"BLEP") = csrmse(m,"%BLEP%","msbias");


* ----------------------------------------------------------------------
*   Analysis of biases
* ----------------------------------------------------------------------

putclose batch "title Analysing biases";
execute "titlebatch";
putclose con / "... Price biases ..." /;

*   Compute descriptive statistics for biases

parameter biasstat(*,a,parvec) Sample mean and variances of biases;

*   ... average bias ...
biasstat("mean",aon,"preg") =
    sum((m,t,i), pmse(m,t,i,aon,"msebias"))
    /(card(m)*card(i)*card(t));

biasstat("mean",aon,"tpc") =
    sum((m,i,j)$AD(i,j), cmse(m,i,j,aon,"msebias"))
    /(card(m)*card(i)*(card(i)-1)/2);

*   ... sum of absolute biases ...
biasstat("SABIAS",aon,"preg") =
    sum((m,t,i), abs(pmse(m,t,i,aon,"msebias")));

biasstat("SABIAS",aon,"tpc")  =
    sum((m,i,j)$AD(i,j), abs(cmse(m,i,j,aon,"msebias")));

*   ... variance of biases ...
biasstat("variance",aon,"preg") =
    sum((m,t,i), sqr(pmse(m,t,i,aon,"msebias")
                    -biasstat("mean",aon,"preg")))
    /(card(m)*card(i)*card(t)-1);

biasstat("variance",aon,"tpc") =
    sum((m,i,j)$AD(i,j), sqr(cmse(m,i,j,aon,"msebias")
                            -biasstat("mean",aon,"tpc")))
    /(card(m)*card(i)*(card(i)-1)/2-1);

putclose con / "... computing medians ..." /;

*   ... median bias over k times m observations for prices (use gdx-rank)

count = 0;
loop(obs$sameas(obs,"o1"),
    loop((m,t,i),
        someSort1(obs+count) = pmse(m,t,i,"%TRAD%","msebias");
        someSort2(obs+count) = pmse(m,t,i,"%BLEP%","msebias");
        count=count+1;
    );
);
uobs(obs) = yes $ (ord(obs) lt count);

pct("median") = 50.0;
$libinclude rank.gms someSort1 uobs someRank pct
biasstat("median","%TRAD%","preg") = pct("median");
pct("median") = 50.0;
$libinclude rank.gms someSort2 uobs someRank pct
biasstat("median","%BLEP%","preg")  =pct("median");

*   ... median bias over k times m observations for costs (use gdx-rank)

count = 0;
loop(obs$sameas(obs,"o1"),
    loop((m,i,j)$AD(i,j),
        someSort1(obs+count) = cmse(m,i,j,"%TRAD%","msebias");
        someSort2(obs+count) = cmse(m,i,j,"%BLEP%","msebias");
        count=count+1;
    );
);
uobs(obs) = yes $ (ord(obs) lt count);

pct("median") = 50.0;
$libinclude rank.gms someSort1 uobs someRank pct
biasstat("median","%TRAD%","tpc") = pct("median");
pct("median") = 50.0;
$libinclude rank.gms someSort2 uobs someRank pct
biasstat("median","%BLEP%","tpc")  =pct("median");

display biasstat;


*   Summarize t-statistics
putclose con / "... computing t-tests for biases ..." /;

*   Double sided test if average price differs from zero, based on
*   asymptotic normality of averages (Lindberg-Levy Central Limit Theorem)

tsummary(aon,alpha) =
    sum((m,t,i) $ [2*(1-errorf(abs[ptval(m,t,i,aon)])) le pcrit(alpha)],1);
tsummary(aon,'neg') =
    sum((m,t,i) $ (ptval(m,t,i,aon) lt 0),1);

display "Number of rejections of null hypothesis for price estimates:"
        ,tsummary;

*   Test for average of cost estimates equals true cost

tsummary(aon,alpha) =
    sum((m,i,j)$ (AD(i,j)
             and [2*(1-errorf(abs[ctval(m,i,j,aon)])) le pcrit(alpha)]),1);
tsummary(aon,"neg") =
    sum((m,i,j)  $ (AD(i,j) and (ctval(m,i,j,aon) lt 0)),1);

display "Number of rejections of null hypothesis for cost estimates:"
        , tsummary;

option kill = tsummary;


* ----------------------------------------------------------------------
*   Analysis of variances
* ----------------------------------------------------------------------

putclose batch "title Analysing variances";
execute "titlebatch";
putclose con / "... pooling variances ..." /;

parameter msestat(totmse,*,aon,parvec);

msestat(totmse,"mean",aon,"preg") = sum(m, psrmse(m,aon,totmse)) / card(m);
msestat(totmse,"mean",aon,"tpc")  = sum(m, csrmse(m,aon,totmse)) / card(m);

msestat(totmse,"variance",aon,"preg") =
    sum(m, sqr(psrmse(m,aon,totmse)-msestat(totmse,"mean",aon,"preg")))
    / (card(m)-1);
msestat(totmse,"variance",aon,"tpc")  =
    sum(m, sqr(csrmse(m,aon,totmse)-msestat(totmse,"mean",aon,"tpc")))
    / (card(m)-1);

display msestat;


* ----------------------------------------------------------------------
*   Tests: 1) Variances of prices are heterogeneous for TRAD and BLEP
*          2) Variances of costs are heterogeneous for BLEP but not TRAD
*          3) BLEP is MSE-dominant over TRAD
*       - Plot of variances
*       - Bartlett test
*       - Levene test (if normality does not hold)
* ----------------------------------------------------------------------

putclose con / "... making vectors suitable for plotting ..." /;

*   Put variances on parameter suitable for plotting

count = 0;
loop(obs$sameas(obs,"o1"),
    loop((m,t,i),
        plotphomo(obs+count,"TRAD") = pmse(m,t,i,"%TRAD%","msevar");
        plotphomo(obs+count,"BLEP") = pmse(m,t,i,"%BLEP%","msevar");
        count=count+1;
    );
);

count = 0;
loop(obs$sameas(obs,"o1"),
    loop((m,i,j)$AD(i,j),
        plotchomo(obs+count,"TRAD") = cmse(m,i,j,"%TRAD%","msevar");
        plotchomo(obs+count,"BLEP") = cmse(m,i,j,"%BLEP%","msevar");
        count=count+1;
    );
);

putclose batch "title Computing Bartlett statistics";
execute "titlebatch";
putclose con / "... testing for equality of variances (Bartlett) ..." /;

bartlett(m,aon,"preg") =
    ((card(i)*card(t)*(card(n)-1))*log(psrmse(m,aon,"mvar"))
     -sum((t,i), (card(n) - 1)*log(pmse(m,t,i,aon,"msevar"))))
    /(1 + (1/(3*(card(i)*card(t)-1)))
        * (sum((t,i), 1/(card(n) - 1))
                    - 1/(card(i)*card(t)*(card(n)-1))));

bartlett(m,aon,"tpc")  =
    ((card(i)*(card(i)-1)/2*(card(n)-1))*log(csrmse(m,aon,"mvar"))
      -sum((i,j)$AD(i,j), (card(n) - 1)*log(cmse(m,i,j,aon,"msevar"))))
    /(1 + (1/(3*(card(i)*(card(i)-1)/2-1)))
        * (sum((i,j)$AD(i,j), 1/(card(n) - 1))
                            - 1/(card(i)*(card(i)-1)/2*(card(n)-1))));

parameter plevYi(m,t,i,a)    The Yi hat parameter in the levene statistic;
parameter plevZij(m,n,t,i,a) The Zij parameter in the levene statistic;
parameter plevZi(m,t,i,a)    The Zi hat parameter in the levene statistic;
parameter plevZ(m,a)         The Z hat parameter in the levene statistic;

parameter clevYi(m,i,j,a)    The Yi hat parameter in the levene statistic;
parameter clevZij(m,n,i,j,a) The Zij parameter in the levene statistic;
parameter clevZi(m,i,j,a)    The Zi hat parameter in the levene statistic;
parameter clevZ(m,a)         The Z hat parameter in the levene statistic;

*   Mean (1) or median (0) for levene measure? Median should be better for
*   skew distributions. Trial shows no difference.

scalar levenemean "Use mean (1) or median (0)" /0/;
plevYi(m,t,i,aon) $ levenemean    = pmean(m,t,i,aon);
clevYi(m,i,j,aon) $ levenemean    = cmean(m,i,j,aon);

parameter levSort(n);
parameter levRank(n);


if((not levenemean),
    putclose con / "... computing medians ..." /
);


loop((m,aon) $ (not levenemean),
    putclose batch "title Computing median price: ",m.tl;
    execute "titlebatch";
    putclose con / "... medians for prices in model ", m.tl /;
    loop((t,i),
        pct("median") = 50.0;
        levSort(n) = pest(m,n,t,i,aon);
$libinclude rank.gms levSort n levRank pct
        plevYi(m,t,i,aon) = pct("median");
    );
    putclose con / "... medians for transport costs in model ", m.tl /;
    putclose batch "title Computing median cost: ",m.tl;
    execute "titlebatch";
    loop((i,j)$AD(i,j),
        pct("median") = 50.0;
        levSort(n) = cest(m,n,i,j,aon);
$libinclude rank.gms levSort n levRank pct
        clevYi(m,i,j,aon) = pct("median");
    );
);

putclose con / "... free some memory ..." /;
solve MDUM using NLP maximising VDUM;

putclose batch "title Computing Levene statistics";
execute "titlebatch";

putclose con / "... computing Levene statistics for prices ..." /;

plevZij(m,n,t,i,aon)  = abs(pest(m,n,t,i,aon)-plevYi(m,t,i,aon));
plevZi(m,t,i,aon)     = sum(n, plevZij(m,n,t,i,aon)) / nn;
plevZ(m,aon)          = sum((t,i), plevZi(m,t,i,aon)) / np;

putclose con / "... free some memory ..." /;
option kill = pest;
solve MDUM using NLP maximising VDUM;
putclose con / "... computing Levene statistics for costs ..." /;

clevZij(m,n,i,j,aon)  = abs(cest(m,n,i,j,aon)-clevYi(m,i,j,aon));
clevZi(m,i,j,aon)     = sum(n, clevZij(m,n,i,j,aon)) / nn;
clevZ(m,aon)          = sum((i,j)$AD(i,j), clevZi(m,i,j,aon)) / nc;

putclose con / "... free some memory ..." /;
option kill = cest;
solve MDUM using NLP maximising VDUM;
putclose con / "... finalize Levene statistics ..." /;

levene(m,aon,"preg") =
    ((np*(card(n)-1))
     *sum((t,i), card(n)*sqr(plevZi(m,t,i,aon)-plevZ(m,aon))))
    /((np-1)*sum((t,i,n), sqr(plevZij(m,n,t,i,aon)-plevZi(m,t,i,aon))));

levene(m,aon,"tpc")  =
    ((nc*(card(n)-1))
     *sum((i,j)$AD(i,j), card(n)*sqr(clevZi(m,i,j,aon)-clevZ(m,aon))))
    /((nc-1)*sum((i,j,n)$AD(i,j), sqr(clevZij(m,n,i,j,aon)
                                     -clevZi(m,i,j,aon))));

display bartlett,levene;

set homtest /bartlett,levene/;
parameter heterostats(*,alpha,parvec,aon) "Summary of Levene and Bartlett";

table critval(parvec,alpha,homtest)
                bartlett    levene
preg.alpha950   16.92       1.88
preg.alpha990   21.67       2.43
tpc.alpha950    61.66       1.39
tpc.alpha990    69.96       1.61    ;


heterostats("bartlett",alpha,parvec,aon) $ (not sameas(parvec,"stc")) =
   sum(m $ (bartlett(m,aon,parvec) ge critval(parvec,alpha,"bartlett")),1);

heterostats("levene",alpha,parvec,aon)  $ (not sameas(parvec,"stc"))  =
   sum(m $ (levene(m,aon,parvec) ge critval(parvec,alpha,"levene")),1);

display heterostats;

*$exit
* ----------------------------------------------------------------------
*   PostScript output for paper
* ----------------------------------------------------------------------

*   Make a latex table summarising the price bias info

file pribias /'%repdir%\tabpricebias.tex'/;

pribias.nd=3;
put pribias;
put "\begin{tabular}{l|c|c|}"             /;
put "          &   TRAD    &   BLEP   \\" /;
put "\hline"                              /;
put "mean      &", biasstat("mean","%TRAD%","preg"),     "  & "
                 , biasstat("mean","%BLEP%","preg"),    " \\ " /;
put "variance  &", biasstat("variance","%TRAD%","preg"), "  & "
                 , biasstat("variance","%BLEP%","preg")," \\ " /;
put "median    &", biasstat("median","%TRAD%","preg"),   "  & "
                 , biasstat("median","%BLEP%","preg"),  " \\ " /;
put "SABIAS    &", biasstat("SABIAS","%TRAD%","preg"),   "  & "
                 , biasstat("SABIAS","%BLEP%","preg"),  " \\ " /;
put "\hline"        /;
put "\end{tabular}" /;
putclose pribias;

*   Make a latex table summarising the cost bias info

file tpcbias /'%repdir%\tabcostbias.tex'/;

tpcbias.nd=3;
put tpcbias;
put "\begin{tabular}{l|c|c|}"             /;
put "          &   TRAD    &   BLEP   \\" /;
put "\hline"                              /;
put "mean      &", biasstat("mean","%TRAD%","tpc"),     "  & "
                 , biasstat("mean","%BLEP%","tpc"),    " \\ " /;
put "variance  &", biasstat("variance","%TRAD%","tpc"), "  & "
                 , biasstat("variance","%BLEP%","tpc")," \\ " /;
put "median    &", biasstat("median","%TRAD%","tpc"),   "  & "
                 , biasstat("median","%BLEP%","tpc"),  " \\ " /;
put "SABIAS    &", biasstat("SABIAS","%TRAD%","tpc"),   "  & "
                 , biasstat("SABIAS","%BLEP%","tpc"),  " \\ " /;
put "\hline"        /;
put "\end{tabular}" /;
putclose tpcbias;


$setglobal GPSTYLE

$setglobal gp_opt0  set data style points
$setglobal gp_opt1  set terminal postscript portrait 'Times-Italic' 10
$setglobal gp_opt2  set size 0.65,0.25
$setglobal gp_opt3  set pointsize 0.5
$setglobal gp_opt4  set key outside
$setglobal gp_opt5  set key width 1
$setglobal gp_opt6  set key samplen 0
$setglobal gp_opt11 set title
$setglobal batch yes


$setglobal gp_opt7  set output '%repdir%\pmsbias.eps'
$setglobal gp_opt8  set xlabel 'Model number'
$setglobal gp_opt9  set ylabel 'Mean squared bias'
$setglobal gp_opt10 set yrange [0:*]

$if %PAPEROUT%==ON $libinclude plot plotpbias m


$setglobal gp_opt7  set output '%repdir%\pmvar.eps'
$setglobal gp_opt8  set xlabel 'Model number'
$setglobal gp_opt9  set ylabel 'Mean variance'
$setglobal gp_opt10 set yrange [0:*]

$if %PAPEROUT%==ON $libinclude plot plotpvar m


$setglobal gp_opt7  set output '%repdir%\cmsbias.eps'
$setglobal gp_opt8  set xlabel 'Model number'
$setglobal gp_opt9  set ylabel 'Mean squared bias'
$setglobal gp_opt10 set yrange [*:*]

$if %PAPEROUT%==ON $libinclude plot plotcbias m


$setglobal gp_opt7  set output '%repdir%\cmvar.eps'
$setglobal gp_opt8  set xlabel 'Model number'
$setglobal gp_opt9  set ylabel 'Mean variance'
$setglobal gp_opt10 set yrange [*:*]

$if %PAPEROUT%==ON $libinclude plot plotcvar m


$setglobal gp_opt7  set output '%repdir%\pmmse.eps'
$setglobal gp_opt8  set xlabel 'Model number'
$setglobal gp_opt9  set ylabel 'Mean MSE'

$if %PAPEROUT%==ON $libinclude plot plotpmse m


$setglobal gp_opt7  set output '%repdir%\cmmse.eps'
$setglobal gp_opt8  set xlabel 'Model number'
$setglobal gp_opt9  set ylabel 'Mean MSE'

$if %PAPEROUT%==ON $libinclude plot plotcmse


    uobs(obs) = yes$(ord(obs) le 200);

$setglobal gp_opt3  set pointsize 0.4
$setglobal gp_opt7  set output '%repdir%\pricehom.eps'
$setglobal gp_opt8  set xlabel 'Estimated regional prices'
$setglobal gp_opt9  set ylabel 'Variance'

$if %PAPEROUT%==ON $libinclude plot plotphomo uobs


$setglobal gp_opt7  set output '%repdir%\costhom.eps'
$setglobal gp_opt8  set xlabel 'Estimated transport cost items'
$setglobal gp_opt9  set ylabel 'Variance'

$if %PAPEROUT%==ON $libinclude plot plotchomo uobs



$label poster
$if %POSTEROUT%==OFF $goto finito
* ----------------------------------------------------------------------
*   PostScript output for PARMA poster
* ----------------------------------------------------------------------

    putclose con / "... PostScript output for Parma poster ..." /;

$setglobal gp_opt1  set terminal postscript portrait 'Helvetica-Oblique' 28
$setglobal gp_opt2  set size 1.1,0.80
$setglobal gp_opt3  set pointsize 1.5
$setglobal gp_opt4  unset key
$setglobal gp_opt12 set lmargin 5
$setglobal gp_opt16 set terminal postscript color color

* SRMSE for prices

$setglobal gp_opt7  set output '%repdir%\plotpmse.eps'
$setglobal gp_opt8  set xlabel 'Model number'
$setglobal gp_opt9  unset ylabel
$setglobal gp_opt11  set title 'Mean MSE of price estimates'
*$setglobal gp_opt11 set title 'Sum of Root Mean Square Errors (SRMSE) for price estimates'
*$setglobal gp_opt11 set title

$libinclude plot plotpmse m

* SRMSE for costs

$setglobal gp_opt7  set output '%repdir%\plotcmse.eps'
$setglobal gp_opt8  set xlabel 'Model number'
$setglobal gp_opt11 set title 'Mean MSE of cost estimates'
*$setglobal gp_opt11 set title 'Sum of Root Mean Square Errors (SRMSE) for cost estimates'

$libinclude plot plotcmse m

* Variance of prices

$setglobal gp_opt7  set output '%repdir%\plotpvar.eps'
$setglobal gp_opt8  set xlabel 'Model number'
$setglobal gp_opt11 set title 'Mean variance of price estimates'
*$setglobal gp_opt11 set title

$libinclude plot plotpvar m

* Variance of costs

$setglobal gp_opt7  set output '%repdir%\plotcvar.eps'
$setglobal gp_opt8  set xlabel 'Model number'
$setglobal gp_opt11 set title 'Mean variance of cost estimates'
*$setglobal gp_opt11 set title

$libinclude plot plotcvar m

* Bias of prices

$setglobal gp_opt4  set key right
$setglobal gp_opt13 set key box
$setglobal gp_opt7  set output '%repdir%\plotpbias.eps'
$setglobal gp_opt8  set xlabel 'Model number'
$setglobal gp_opt11 set title 'MSBIAS for price estimates'
$setglobal gp_opt10 set yrange [*:*]
*$setglobal gp_opt11 set title

$libinclude plot plotpbias m

* Bias of costs

$setglobal gp_opt4  unset key
$setglobal gp_opt7  set output '%repdir%\plotcbias.eps'
$setglobal gp_opt8  set xlabel 'Model number'
$setglobal gp_opt11 set title 'MSBIAS for transport cost estimates'
$setglobal gp_opt10 set yrange [*:*]
*$setglobal gp_opt11 set title

$libinclude plot plotcbias m

* Hetero(homo)geneity of prices

$setglobal gp_opt3  set pointsize 0.5
$setglobal gp_opt7  set output '%repdir%\plotphomo.eps'
$setglobal gp_opt8  set xlabel 'Estimated prices'
$setglobal gp_opt11 set title 'Sample variance of price estimates'
$setglobal gp_opt15 set xtics 250

    uobs(obs) = yes$(ord(obs) le card(m)*card(i));

* Dangerous!!!! Watch out! Series are interchanged to make plot nicer! Turned right with java program afterwards.
    parameter plottemp(obs);

    plottemp(obs)         = plotphomo(obs,"BLEP");
    plotphomo(obs,"BLEP") = plotphomo(obs,"TRAD");
    plotphomo(obs,"TRAD") = plottemp(obs);

    plottemp(obs)         = plotchomo(obs,"BLEP");
    plotchomo(obs,"BLEP") = plotchomo(obs,"TRAD");
    plotchomo(obs,"TRAD") = plottemp(obs);
    plottemp(obs) = 0;

$libinclude plot plotphomo uobs

* Hetero(homo)geneity of costs

$setglobal gp_opt7  set output '%repdir%\plotchomo.eps'
$setglobal gp_opt8  set xlabel 'Estimated transport costs'
$setglobal gp_opt11 set title 'Sample variance of cost estimates'
$setglobal gp_opt15 set xtics 1000


    uobs(obs) = yes$(ord(obs) le card(m)*card(i)*(card(i)-1)/2);
*    uobs(obs) = yes$(ord(obs) le 2250);
$libinclude plot plotchomo uobs

file repl /repl.bat/;

put repl 'REM Replace colours in EPS files';

put / 'cd %repdir%'/;
put / 'java repStr plotphomo.eps " Crs" " blaha"';
put / 'java repStr plotphomo.eps " Pls" " Crs"';
put / 'java repStr plotphomo.eps " blaha" " Pls"';

put / 'java repStr plotchomo.eps " Crs" " blaha"';
put / 'java repStr plotchomo.eps " Pls" " Crs"';
put / 'java repStr plotchomo.eps " blaha" " Pls"';


put / 'java repStr plotpmse.eps "/LT0 { PL [] 1 0 0 DL } def" "/LT0 { PL [] 0.6 0 0 DL } def"';
put / 'java repStr plotcmse.eps "/LT0 { PL [] 1 0 0 DL } def" "/LT0 { PL [] 0.6 0 0 DL } def"';
put / 'java repStr plotpbias.eps  "/LT0 { PL [] 1 0 0 DL } def" "/LT0 { PL [] 0.6 0 0 DL } def"';
put / 'java repStr plotcbias.eps  "/LT0 { PL [] 1 0 0 DL } def" "/LT0 { PL [] 0.6 0 0 DL } def"';
put / 'java repStr plotpvar.eps   "/LT0 { PL [] 1 0 0 DL } def" "/LT0 { PL [] 0.6 0 0 DL } def"';
put / 'java repStr plotcvar.eps   "/LT0 { PL [] 1 0 0 DL } def" "/LT0 { PL [] 0.6 0 0 DL } def"';
put / 'java repStr plotphomo.eps  "/LT0 { PL [] 1 0 0 DL } def" "/LT0 { PL [] 0 0.6 0 DL } def"';
put / 'java repStr plotchomo.eps  "/LT0 { PL [] 1 0 0 DL } def" "/LT0 { PL [] 0 0.6 0 DL } def"';

put / 'java repStr plotpmse.eps "/LT1 { PL [4 dl 2 dl] 0 1 0 DL } def" "/LT1 { PL [4 dl 2 dl] 0 0.6 0 DL } def"';
put / 'java repStr plotcmse.eps "/LT1 { PL [4 dl 2 dl] 0 1 0 DL } def" "/LT1 { PL [4 dl 2 dl] 0 0.6 0 DL } def"';
put / 'java repStr plotpbias.eps  "/LT1 { PL [4 dl 2 dl] 0 1 0 DL } def" "/LT1 { PL [4 dl 2 dl] 0 0.6 0 DL } def"';
put / 'java repStr plotcbias.eps  "/LT1 { PL [4 dl 2 dl] 0 1 0 DL } def" "/LT1 { PL [4 dl 2 dl] 0 0.6 0 DL } def"';
put / 'java repStr plotpvar.eps   "/LT1 { PL [4 dl 2 dl] 0 1 0 DL } def" "/LT1 { PL [4 dl 2 dl] 0 0.6 0 DL } def"';
put / 'java repStr plotcvar.eps   "/LT1 { PL [4 dl 2 dl] 0 1 0 DL } def" "/LT1 { PL [4 dl 2 dl] 0 0.6 0 DL } def"';
put / 'java repStr plotphomo.eps  "/LT1 { PL [4 dl 2 dl] 0 1 0 DL } def" "/LT1 { PL [4 dl 2 dl] 0.6 0 0 DL } def"';
put / 'java repStr plotchomo.eps  "/LT1 { PL [4 dl 2 dl] 0 1 0 DL } def" "/LT1 { PL [4 dl 2 dl] 0.6 0 0 DL } def"';

putclose / ;

execute 'repl.bat';

$label finito