RASing a matrix
Posted: Thu May 28, 2020 10:22 am
How do I RASe a matrix in GAMS?
Code: Select 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;