$title Simple investment example with varying weight for risk using GUSS (GUSSRISK,SEQ=399)

$onText
Portfolio selection model solved for alternative risk aversion parameters.

Keywords: nonlinear programming, portfolio optimization, GUSS, scenario analysis,
          investment planning
$offText

Set
   stocks 'potential investments'                  / buystock1*buystock4 /
   events 'equally likely return states of nature' / event1*event10      /;

Alias (stocks,s,sp), (events,e,ep);

Parameter
   prices(stocks) 'purchase prices of the stocks' / buystock1  22
                                                    buystock2  30
                                                    buystock3  28
                                                    buystock4  26 /
   funds 'total investable funds' / 500 /;

Table returns(events,stocks) 'returns by state of nature event'
           buystock1  buystock2  buystock3  buystock4
   event1          7          6          8          5
   event2          8          4         16          6
   event3          4          8         14          6
   event4          5          9         -2          7
   event5          6          7         13          6
   event6          3         10         11          5
   event7          2         12         -2          6
   event8          5          4         18          6
   event9          4          7         12          5
   event10         3          9         -5          6;

Parameter
   mean (stocks)        'mean returns to x(stocks)'
   covar(stocks,stocks) 'variance covariance matrix';

mean(s)     = sum(e, returns(e,s))/card(e);
covar(s,sp) = sum(e,(returns(e,sp) - mean(sp))*(returns(e,s) - mean(s)))/card(e);

Scalar rap 'risk aversion parameter' / 0 /;

Variable
   invest(stocks) 'money invested in each stock'
   obj            'number to be maximized';

Positive Variable invest;

Equation
   objj     'objective function'
   investav 'investment funds available';

objj.. obj =e= sum(s, mean(s)*invest(s)) - rap*(sum((s,sp), invest(s)*covar(s,sp)*invest(sp)));

investav.. sum(s, prices(s)*invest(s)) =l= funds;

Model evportfol / all /;

solve evportfol using nlp maximizing obj;

Scalar variance;
variance = sum(s, sum(sp, invest.l(s)*covar(s,sp)*invest.l(sp)));

display variance;

Set rapscenarios 'risk aversion parameters' / r0*r25 /;

Parameter
   riskaver(rapscenarios) 'risk aversion coeficient by risk aversion parameter'
                          / r0  0.00000, r1  0.00025, r2  0.00050, r3  0.00075
                            r4  0.00100, r5  0.00150, r6  0.00200, r7  0.00300
                            r8  0.00500, r9  0.01000, r10 0.01100, r11 0.01250
                            r12 0.01500, r13 0.02500, r14 0.05000, r15 0.10000
                            r16 0.30000, r17 0.50000, r18 1.00000, r19 2.50000
                            r20 5.00000, r21 10.0000, r22 15.    , r23 20.
                            r24 40.    , r25 80.                              /;

Parameter
   stockoutput(rapscenarios,stocks) 'results for invest based on model runs with varying rap'
   objlevel(rapscenarios)           'results for objective function from model runs with varying rap'
   investavshadow(rapscenarios)     'results for funds shadow price from model runs with varying rap';

Set dict / rapscenarios.scenario. ''
           rap         .param   .riskaver
           invest      .level   .stockoutput
           obj         .level   .objlevel
           investav    .marginal.investavshadow /;

solve evportfol using nlp maximizing obj scenario dict;

display stockoutput;

Parameter output(*,rapscenarios);

loop(rapscenarios,
   rap = riskaver(rapscenarios);
*  load in solution information
   invest.l(s) = stockoutput(rapscenarios,s);
   obj.l       = objlevel(rapscenarios);
   investav.m  = investavshadow(rapscenarios);
*  compute some items
   investav.l = sum(s,invest.l(s));
   variance   = sum((s,sp), invest.l(s)*covar(s,sp)*invest.l(sp));
   output('rap'      ,rapscenarios) = rap;
   output(stocks     ,rapscenarios) = invest.l(stocks);
   output('obj'      ,rapscenarios) = obj.l;
   output('mean'     ,rapscenarios) = sum(s, mean(s)*invest.l(s));
   output('var'      ,rapscenarios) = variance;
   output('std'      ,rapscenarios) = sqrt(variance);
   output('shadprice',rapscenarios) = investav.m;
   output('idle'     ,rapscenarios) = funds - investav.l;
);
display output;