(!******************************************************
Mosel Example Problems
======================
file 021folio_robust.mos
````````````````````````
Single period portfolio allocation problem.
The present model produces a selection of
assets that is robust against realization of
worst case with a high probability.
Features:
- Using the 'ellipsoid' uncertainty set
- Retrieving the worst value of an uncertain variable
What to look at:
- Selection of highest worst case value asset,
but low expected value for the problem
with highest protection level
- Selection of highest exepected value asset,
but low worst case value for the problem
with high protection level
- Selection of a diversied set of assets for
various value of the protection level.
- Using worst value of an uncertain variable as
an indicator for reoptimizing the portfolio
Model is inspired by work from the following book:
A. Ben-Tal, L. El Ghaoui, and A. Nemirovski. Robust Optimization.
Princeton Series in Applied Mathematics, 2009.
(c) 2014 Fair Isaac Corporation
author: S. Lannez, Mar. 2014, rev. Dec. 2014
*******************************************************!)
model "Robust portfolio optimization"
uses "mmrobust", "random"
uses "mmsystem"
parameters
ZEROTOL=1e-6
SEED=12345
NMC=500000
N=1.5 ! Worst case metric
end-parameters
declarations
Problems: set of mpproblem
ProtectLevel: set of real
mp_problemA: mpproblem ! Worst case optimization
mp_problemB: array(ProtectLevel) of mpproblem ! Robust optimization
NSHARES = 25 ! Max number of shares
Shares = 1..NSHARES ! Set of shares
PRICE: array(Shares) of real ! Expected return on investment
VAR: array(Shares) of real ! Uncertainty measure (deviation) per SHARE
expReturn: mpvar ! Expected portfolio value
wstReturn: mpvar ! Worst case portfolio value
frac: array(Shares) of mpvar ! Fraction of capital used per share
frac_sol: array(Shares,Problems) of real
obj: mpvar
e: array(Shares) of uncertain ! Deviation of share values
PRICEthr: array(Problems,Shares) of real ! The price threshold
end-declarations
!***************************Configuration*************************
setparam("XPRS_LOADNAMES",true)
!***************************Subroutines***************************
!**** Create the nominal model ****
procedure create_nominalmodel
! Nominal model
expReturn = (sum(s in Shares) PRICE(s)*frac(s))
wstReturn = sum(s in Shares) ( PRICE(s) - N*VAR(s) )*frac(s)
! Spend all the budget
sum(s in Shares) frac(s) = 1
end-procedure
!**** Optimize for the worst case realization ****
procedure solve_det
obj = wstReturn
maximize(obj)
end-procedure
!**** Optimize for a given protection level ****
procedure solve_rob(k: real)
! The value variation domain
assert(k>=0 and k<=1)
sum(s in Shares) (e(s) / (N*VAR(s)))^2 <= (k)^2
! The robust constraint
obj <= sum(s in Shares) (PRICE(s) + e(s)) * frac(s)
maximize(obj)
end-procedure
!***************************Main model***************************
!**** Input data generation ****
setmtrandseed(12345)
cnt := 0
forall(s in Shares, cnt as counter) do
PRICE(s) := round((100*(NSHARES-cnt+1)/NSHARES))
c := (1+sqrt((NSHARES-s)/NSHARES))/3
VAR(s) := round(PRICE(s)*c*100)/100
end-do
writeln("\n *** Input Data *** ")
writeln
writeln("Shares | Mean | S.Dev. | Worst value")
forall(s in Shares)
writeln(strfmt(s,6), " |", strfmt(PRICE(s),5), " |", strfmt(VAR(s),7,2),
" |", strfmt(PRICE(s)-N*VAR(s),6,3))
writeln
!**** Optimize the worst case ****
with mp_problemA do
create_nominalmodel
solve_det
forall(s in Shares) frac_sol(s,mp_problemA) := frac(s).sol*100
! Price set by opponent is worst case price
forall(s in Shares) PRICEthr(mp_problemA,s) := PRICE(s) - N*VAR(s)
end-do
!**** Optimize the 'budgeted' worst case ****
Ks := [0, ! No variation on average
0.01, ! +/- 1% variation on the list
0.10, ! +/- 10% variation on the list
0.25, ! +/- 25%
1] ! +/-100% variation on the list
forall(k in Ks) do
create(mp_problemB(k))
with mp_problemB(k) do
create_nominalmodel
solve_rob(k)
forall(s in Shares) frac_sol(s,mp_problemB(k)) := frac(s).sol*100
forall(s in Shares) PRICEthr(mp_problemB(k),s) := (PRICE(s) + getsol(e(s)))
end-do
end-do
!**** Print results ****
writeln("\n *** Portfolio allocation results *** ")
write("\n | protection level")
write("\nShares | Wst price | Price | A |")
forall(k in Ks) write(" ", strfmt(k*100,3), "% |" ) ; writeln
forall(s in Shares) do
write(strfmt(s,6), " |", strfmt(PRICE(s)-N*VAR(s),10,0), " |",
strfmt(PRICE(s),6), " | ")
forall(mp in Problems) do
if (frac_sol(s,mp)>1) then
write(strfmt(frac_sol(s,mp),3,0), "% | ")
else
write(" | ")
end-if
end-do
writeln
end-do
writeln("\n *** Worst case price *** ")
write("\n | protection level")
write("\nShares | Wst price | Price | A |")
forall(k in Ks) write(" ", strfmt(k*100,3), "% |" ) ; writeln
forall(s in Shares) do
write(strfmt(s,6), " |", strfmt(PRICE(s)-N*VAR(s),10,0), " |",
strfmt(PRICE(s),6), " | ")
forall(mp in Problems) do
if (frac_sol(s,mp)>1) then
write(strfmt(PRICEthr(mp,s),3,0), "* | ")
else
write(strfmt(PRICEthr(mp,s),4,0), " | ")
! write(" | ")
end-if
end-do
writeln
end-do
!**** Simulate results (Monte-Carlo simulation) ****
writeln("\nRunning Monte-Carlo simulation...")
Values := {0,20,40,60,80,100}
declarations
cntV: dynamic array(Problems,Values) of real
s1: real
s2: real
end-declarations
forall(mp in Problems) do
expRev(mp) := 0.0 ; wstRev(mp) := 0.0 ; expDev(mp) := 0.0
forall(v in Values) cntV(mp,v) := 0
c := 0.0 ! Number of draws
s1 := 0 ; s2 := 0;
forall(i in 1..NMC, c as counter) do
totalVal := 0.0 ; totalRisk := 0.0
forall(s in Shares | frac_sol(s,mp)>0) do
! draw a price
p := normal(PRICE(s),VAR(s))
! calculate share revenue
r := frac_sol(s,mp)*p/100
! calculate total revenue
totalVal += r ! summing up revenue
end-do
forall(v in Values) do
if (totalVal>v) then
cntV(mp,v) += 1
end-if
end-do
expRev(mp) += totalVal
s1 += totalVal
s2 += totalVal^2
end-do
expRev(mp) := expRev(mp) / c
expDev(mp) := sqrt(NMC*s2 - s1^2)/(NMC)
forall(v in Values) cntV(mp,v) := cntV(mp,v) / c
end-do
!**** Print simulation results ****
write("\n | protection level")
write("\n | A |")
forall(k in Ks) write(" ", strfmt(k*100,3), "% |" )
write("\n Expected value |")
forall(mp in Problems) write(strfmt(round(expRev(mp)),5), " |")
write("\n Standard deviation |")
forall(mp in Problems) write(strfmt(round(expDev(mp)),5), " |")
writeln
forall(v in Values) do
write("\n P (value>"+textfmt(v,3)+" |")
forall(mp in Problems) write(strfmt(cntV(mp,v),5,2) , " |")
end-do
end-model
|
