Implementation

The following Mosel model implements the mathematical model from the previous section. Notice how the ellipsoid uncertainty set is added to the original problem to create a robust optimization problem.

Also, take a look at how the Monte-Carlo simulation method is applied to quantify the quality of the solution: in NMC=5000 iterations we draw random values for the expected value of every asset by applying a normal distribution centered around its price (mean value) and using its known variance. We count the occurrence of objective function values (or more precisely, which value range the resulting solution value belongs to) to determine the probabilities of different solution qualities. This method is particularly helpful to gain some insight about the solution quality when it is difficult to determine the exact shape of the distribution function of a random variable.

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