User Tools

Site Tools


gams:rasing_a_matrix

RASing a Matrix

A frequent request is for code samples for RASing a matrix. Below is a standalone program that (in this case) RASes a consumption flows matrix. With appropriate modifications, it can be used for square matrices (IO tables or SAMs) as well. This has been contributed by Tom Rutherfold and he uses it frequently for matrices up to 30×30 without any difficulty at all.

$title Program to carry out RAS if consumption shares

Sets
    i 'sectors'         / 'Farm Food Crops'
                          'Farm Nonfood Crops'
                          'Livestock'
                          'Forestry'
                          'Fishery'
                          'Coal Gas & Oil'
                          'Other Mining'
                          'Food Bev & Tobacco'
                          'Textiles & Leather'
                          'Wood & Furniture'
                          'Paper & Printing'
                          'Chemicals & Fertilizers'
                          'Refining LNG & Coal'
                          'Nonmetallic Mineral'
                          'Basic Metals'
                          'Metal Prod & Machines'
                          'Other Industry'
                          'Elect Gas & Water'
                          'Construction'
                          'Trade & Storage'
                          'Restaurants & Hotels'
                          'Rail & Road Transport'
                          'Sea Air Transport'
                          'Financial Serv & Insurance'
                          'Real Estate'
                          'Public Admin'
                          'Social & Other Services' /
    hh 'household type' / 'Ag worker'
                          'Ag owners'
                          'Low wage' 
                          'High wage'/;
                          
TABLE rasData(*,*) 'flows and row and col sum controls'
                             'Ag worker' 'Ag owners' 'Low wage' 'High wage'     csum
'Farm Food Crops'                165.56     5356.08    1331.07     2757.31   6394.42
'Farm Nonfood Crops'             149.77      410.60     154.18      357.37    990.61
'Livestock'                      141.53      905.99     501.23     1904.72   2661.24
'Forestry'                        56.70      141.65      55.69       56.52    259.31
'Fishery'                        170.94      849.68     394.32      932.19   1569.93
'Coal Gas & Oil'                    .00         .00        .00         .00       .00
'Other Mining'                      .22         .86        .29         .56       .21
'Food Bev & Tobacco'            1169.48     5245.54    2777.09     6066.52  13911.94
'Textiles & Leather'              82.61      529.25     279.48      806.17   1367.44
'Wood & Furniture'                18.70      125.78      63.12      249.40    329.88
'Paper & Printing'                 5.10       54.76      33.10      151.39    161.59
'Chemicals & Fertilizers'         47.13      432.05     299.22      844.25   1160.43
'Refining LNG & Coal'             56.26      515.80     357.23     1007.92   1385.44
'Nonmetallic Mineral'              2.33       21.38      14.80       41.75     57.41
'Basic Metals'                      .00         .00        .00         .00       .00
'Metal Prod & Machines'           44.89      481.86     291.23     1332.10   1419.86
'Other Industry'                   3.67       39.38      23.80      108.85    116.05
'Elect Gas & Water'                5.63       51.97     119.31      413.36    523.35
'Construction'                      .00         .00        .00         .00       .00
'Trade & Storage'                  3.20       18.64      39.30      113.54   5737.68
'Restaurants & Hotels'           108.44     1045.05    1061.97     2564.35   4553.91
'Rail & Road Transport'           42.85      313.51     600.14     1737.71   3552.82
'Sea Air Transport'               13.25      147.42     277.10      994.59   1607.77
'Financial Serv & Insurance'      23.17      300.03     136.31      589.15    883.65
'Real Estate'                     87.98      779.65     700.31     1979.36   3287.49
'Public Admin'                      .00         .00        .00         .00       .00
'Social & Other Services'        116.65      853.46     987.11     4353.09   9063.46
rsum                            2516.05    18620.37   10497.40    29362.18
;                                                            

Parameter                                                    
    conflow(i,hh) 'Initial private consumption flows'        
    c0(i)         'Control vector'                           
    con(hh)       'aggregate consumption levels';

conflow(i,hh) = rasData(i,hh);
c0(i) = rasData(i,'csum');
con(hh) = rasData('rsum',hh);

alias (i,rr), (hh,cc);

Parameter
    a0(rr,cc)      'Initial coefficients matrix to RAS'
    a1(rr,cc)      'Final coefficients matrix after RAS'
    rasmat0(rr,cc) 'Initial flows matrix to RAS'
    ct(cc)         'RAS column control totals'
    rt(rr)         'RAS row control totals'
    ratio          'Adjustment parameter on control totals'
    checkcol       'Check sum of column control totals'
    checkrow       'Check sum of row control totals'
    sumccc(cc)     'Original column sums of RAS matrix'
    sumrrr(rr)     'Original row sums of RAS matrix';
    
Variables
    DEV            'Deviations'
    RASMAT(rr,cc)  'RASed matrix'
    R1(rr)         'Rho of RAS matrix'
    S1(cc)         'Sigma of RAS matrix'
    LOSS           'Objective (loss) function value';
    
* Parameter initialization
sumccc(cc)     = sum(rr, conflow(rr,cc));
sumrrr(rr)     = sum(cc, conflow(rr,cc));
a0(rr,cc)      = conflow(rr,cc) / sumccc(cc);
rasmat0(rr,cc) = a0(rr,cc) * con(cc);
ct(cc)         = con(cc);
rt(rr)         = c0(rr);
ratio          = sum(rr, rt(rr)) / sum(cc, ct(cc));
ct(cc)         = ct(cc) * ratio;
checkcol       = sum(cc, ct(cc));
checkrow       = sum(rr, rt(rr));
display ratio, checkcol, checkrow;
display conflow, a0;
display con, ct;
display c0, rt;

* Variable initialization
DEV.l = 0.0;
R1.l(rr) = 1;
S1.l(cc) = 1;
RASMAT.l(rr,cc) = a0(rr,cc) * ct(cc);
con(cc) = ct(cc);

Equations
    biprop(rr,cc) 'Bi-proportionality for RAS matrix'
    devsq         'Definition of squared deviations'
    obj           'Objective function'
    rconst(rr)    'Row constraint'
    cconst(cc)    'Column constraint';
    
biprop(rr,cc).. RASMAT(rr,cc) =e= R1(rr)*S1(cc)*rasmat0(rr,cc);
cconst(cc)..    ct(cc) =e= sum(rr, RASMAT(rr,cc));
rconst(rr)..    rt(rr) =e= sum(cc, RASMAT(rr,cc));
devsq..         DEV    =e= sum( (rr,cc)$rasmat0(rr,cc),
                             sqr( (RASMAT(rr,cc) - rasmat0(rr,cc)) / rasmat0(rr,cc)) );
obj..           LOSS   =e= sum(rr, sqr(R1(rr)) + sqr(1/R1(rr)) )
                         + sum(cc, sqr(S1(cc)) + sqr(1/S1(cc)) ) ;

* Variable bounds
RASMAT.lo(rr,cc) = 0;
R1.lo(rr)        = 0.01;
S1.lo(cc)        = 0.01;

Model consumeRAS / all - devsq /;

solve consumeRAS using nlp minimizing LOSS;

display rasmat.l, r1.l, s1.l;
a1(rr,cc) = rasmat.l(rr,cc) / ct(cc);
display a0, a1 ;
IMPRESSUM / LEGAL NOTICEPRIVACY POLICY gams/rasing_a_matrix.txt · Last modified: 2020/05/28 08:26 by Michael Bussieck