How do I RASe a matrix in GAMS?
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;