Valuation of Securities: 
Computational Aspects

Institut für Mathematische Stochastik, Technische Universität Dresden, 
November 5, 2002

Jan Hurt
Division of Financial & Insurance Mathematics
Faculty of Mathematics & Physics
CHARLES UNIVERSITY PRAGUE

Home page: http://adela.karlin.mff.cuni.cz/~hurt/

Introduction

Securities (Assets)

Primary: stocks (property), bonds (debt)

Derivative securities (derivatives): options, forwards, futures, swaps

Vanilla options are simple either European or American options.

Portfolio: a set of assets

N assets 1,  ... , N, the investor's wealth W, usually W = 1, then a portfolio is represented by proportions (Subscript[x, 1]),  ... , (Subscript[x, N]),  (∑_(n = 1)^N) (Subscript[x, n]) = 1. The case (Subscript[x, n]) < 0 represents a short sale.

Return, expected return, risk

(Subscript[P, t]) ... price of an security at time t
Rate of return (return): (Subscript[R, t]) = (((Subscript[P, t + 1]) - (Subscript[P, t]))/(Subscript[P, t]))or (Subscript[r, t]) = ln(
            1 + (Subscript[R, t])) = ln(
              (Subscript[P, t + 1])/(Subscript[P, t])) = :(Subscript[p, t + 1]) - (Subscript[p, t]) (logarithmic prices)
E (Subscript[R, t]) ... the expected return, σ = (Sqrt[var (Subscript[R, t])])... the risk (volatility)

Positions (from the point of view of a trader):
long position - the trader possesses the asset, he/she buys 	      +1
short position - the trader sells the asset to a counterpart    	       -1

Present Value

For a cash flow CF =((Subscript[CF, 1]),  ... , (Subscript[CF, T])) and a  valuation interest rate r, the formula for the present value reads:

PV(CF,r)=(∑_(t = 0)^T) ((Subscript[CF, t])/(((1 + r))^t)).

More generally, CF =((Subscript[CF, (Subscript[t, 1])]),  ... , (Subscript[CF, (Subscript[t, T])])), with the valuation date (Subscript[t, 0]), and the (yearly) interest rates r=((Subscript[r, (Subscript[t, 1])]),  ... , (Subscript[r, (Subscript[t, T])])), often called yield curve

PV(CF, r)=\(TraditionalForm\`∑\+\(i = 1\)\%T CF\_\(t\_i\)\/\((1 + r\_\(t\_i\))\)\^\(t\_i - t\_0\)\).

In fact, this is the scalar product of the cash flow and the vector of discount factors. Implementation heavily uses listability:

cf = {cf0, cf1, cf2, cf3};
len = Length[cf]
disc = 1/(1 + r)^(Range[len] - 1)
cf . disc

pv1[cf_, r_] := 
      cf . (1/(1 + r)^(Range[Length[cf]] - 
                  1))
pv1[cf, r]

Modification for the payments going on given different dates and we evaluate the present value with respect to the startday:

Date[]
today = Date[][[Range[3]]]

Needs["Miscellaneous`Calendar`"]

dif = 
      Map[DaysBetween[#, today] &, {{2004, 3, 28}, {2003, 1, 13}, {2003, 1, 
            1}}]
discounts = 1/((1 + r)^(dif/360))

pv2[cf_, dates_, r_, startday_] := 
      Module[{discfactors, dif}, 
dif = 
          Map[DaysBetween[#, startday] &, 
            dates]; 
discfactors = 
          1/((1 + r)^(dif/360)); 
cf . 
          discfactors]

dates = {{2004, 3, 28}, {2003, 1, 13}, {2003, 1, 
            1}};
cf = {1000, -900, -50};
pv2[cf, dates, 0.1, today]

Yield curves

Use for valuation cash flows, typically bonds

Given the observable (market) interest rates at different times, we need to approximate for times where the interest rates are unknown, not available on the market. Two approaches: a) strict numerical interpolation, b) fitting

Here we have observed market interest rates in the form {days,  yearly interest rate}:

ir3 = {{7, 4.74`}, {14, 4.72833333333333`}, {30, 4.705`}, {60, 
            4.595`}, {90, 4.53`}, {180, 4.42`}, {270, 4.415`}, {360, 
            4.41`}, {730, 4.51`}, {1095, 4.68`}, {1460, 4.85`}, {1825, 
            4.99`}, {2190, 5.09`}, {2555, 5.19`}, {2920, 5.27`}, {3285, 
            5.33`}, {3650, 5.38`}, {4015, 5.42`}, {4380, 5.46`}, {4745, 
            5.48666666666667`}, {5110, 5.51333333333333`}, {5475, 
            5.54`}};

p2ir3 = 
        ListPlot[ir3  ];

Needs["Statistics`Master`"]

Fitting by a polynomial (not recommendable, just for an illustration how badly it works):

(* 
      Clear[(Subscript[β, 0]), (Subscript[β, 1]), (Subscript[β, 2]), (Subscript[β, 3])] *) 
fit5 = 
      NonlinearFit[
        ir3, (Subscript[β, 0]) + (Subscript[β, 1])*t + (Subscript[β, 2])*(t^2) + (Subscript[β, 3])*(t^3), 
        t, {(Subscript[β, 0]), (Subscript[β, 1]), (Subscript[β, 2]), (Subscript[β, 3])}]

p5ir3 = Plot[fit5, {t, 0, 6230}];
Show[p2ir3, p5ir3];

A better approach:

Clear[b0, b1, b2, c]
fit3 = 
      NonlinearFit[ir3, 
        b0 + (b1 + b2)*(1 - Exp[-t/c])/(t/c) - 
          b2*Exp[-t/c], t, {b0, b1, b2, c}]

p1ir3 = Plot[fit3, {t, 0, 6230}];

Show[p1ir3, p2ir3];

Another better approach: Fitting by cubic splines

xxx = (tri3 = 
              Transpose[ir3])[[1]];
uuu = tri3[[2]];
p1 = Table[1, {Length[uuu]}];
plotsplines[xxx, uuu, p1]

Options

Black-Scholes formula and related topics

Simulation of a simple random walk

ListPlot[NestList[(# + Random[] - 0.5) &, 10, 50], 
        PlotJoined -> True];

Next we show 100 simulated simple random walks, each of length 50, simultaneously:

Needs["Statistics`Master`"]
tab1 = 
        Table[ListPlot[NestList[(# + Random[] - 0.5) &, 10, 50], 
            PlotJoined -> True, 
            DisplayFunction -> 
              Identity], {100}];
tab1[[1]];

Show[tab1, DisplayFunction -> $DisplayFunction];

Next we show a simulation of the daily price development of a stock with the expected return 5 per cent p.a. and the volatility 0.003. The daily increments possess a normal distribution.

s0 = 100; μ = 0.05/365; σ = 
      0.003;
distr = 
        NormalDistribution[μ, σ];
tab2 = 
        Table[ListPlot[NestList[#*(1 + Random[distr]) &, s0, 360], 
            PlotJoined -> True, 
            DisplayFunction -> 
              Identity], {30}];
Show[tab2, DisplayFunction -> $DisplayFunction];

Black-Scholes formula

Notation: s ... current price of underlying asset at time t, k ... strike price, r ..  riskless interest rate, T ... expiry date, σ ... volatility of the underlying asset
call (put) ... the value of a CALL (PUT) option

Clear[Φ]
Φ[x_] := 
      (1/2) (1 + Erf[(x/(Sqrt[2]))])
ϕ[x_] := Exp[-x^2/2]/Sqrt[2*π]

d1 = (Log[s/k] + (r + σ^2/2)*(T - t))/(σ*
              Sqrt[T - t]);
d2 = (Log[s/k] + (r - σ^2/2)*(T - t))/(σ*
              Sqrt[T - t]);
call[s_, k_, σ_, T_, t_, r_] := 
      Module[{d1, d2}, 
        d1 = (Log[s/k] + (r + σ^2/2)*(T - t))/(σ*
                Sqrt[T - t]); 
d2 = (Log[
                  s/k] + (r - σ^2/2)*(T - t))/(σ*
                Sqrt[T - t]); 
s*Φ[d1] - 
          k*Exp[-r*(T - t)]*Φ[d2]]
put[s_, k_, σ_, T_, t_, r_] := 
      Module[{d1, d2}, 
        d1 = (Log[s/k] + (r + σ^2/2)*(T - t))/(σ*
                Sqrt[T - t]); 
d2 = (Log[
                  s/k] + (r - σ^2/2)*(T - t))/(σ*
                Sqrt[T - t]); 
k*
            Exp[-r*(T - t)]*Φ[-d2] - 
          s*Φ[-d1]]

Implied volatility: price ... the market price of an option at time t

impliedcall[s_, k_, T_, t_, r_, price_] := 
      sd /. FindRoot[
          call[s, k, sd, T, t, r] == price, {sd, 
            0.2}]
impliedput[s_, k_, T_, t_, r_, price_] := 
      sd /. FindRoot[
          put[s, k, sd, T, t, r] == price, {sd, 0.2}]

Real data: Underlying Deutsche Bank April 23, 1999 - Price EUR 52.44

optionDB46c[pos_, S_] = {pos, 1, 5.26, S, 
          46};
optionDB48c[pos_, S_] = {pos, 1, 3.75, S, 
          48};
optionDB50c[pos_, S_] = {pos, 1, 2.51, S, 
          50};
optionDB46p[pos_, S_] = {pos, -1, 0.74, S, 
          46};
optionDB48p[pos_, S_] = {pos, -1, 1.39, S, 
          48};
optionDB50p[pos_, S_] = {pos, -1, 2.33, S, 50};

call[52.44, 50, 0.2, 20/52, 0, 0.05]
put[49, 50, 0.2, 20/52, 0, 0.05]

Plot[{call[s, 46, 0.3, 28/360, 0, 0.05], 
          call[s, 50, 0.3, 28/360, 0, 0.05]}, {s, 40, 53}, 
        PlotRange -> All];

Plot[{call[  s, 50, 0.2, 22/52, 0, 
            0.05], put[s, 50, 0.2, 20/52, 0, 0.05]}, {s, 0, 70}, 
        PlotRange -> All];

Real data: Volkswagen

{impliedcall[70.72, k1 = 67.5, T = 58/360, 
        t = 0, r = 0.05, c1 = 6.31], 
impliedcall[70.72, 
        k2 = 70, 58/360, 0, 0.05, c2 = 4.92], 
impliedcall[
        70.72, k3 = 72.5, 58/360, 0, 0.05, c3 = 3.37]}

    {impliedput[70.72, 67.5, 58/360, 0, 0.05, 
        p1 = 2.92], 
impliedput[70.72, 70, 58/360, 0, 0.05, 
        p2 = 4.08], 
impliedput[70.72, 72.5, 58/360, 0, 
        0.05, p3 = 5.48]}

Portfolio of options

Butterfly spread, the value of the portfolio at expiry

valuepor[s_] :=   Max[s - k1, 0] + 
        Max[k3 - s, 0] - Max[s - k2, 0] - 
        Max[k2 - s, 0]
Plot[valuepor[s], {s, 40, 100}, PlotRange -> All];

Butterfly spread valued by BS formula:

Plot[
        call[  s, k1, 0.2, 22/52, 0, 0.05] - 
          call[s, k2, 0.2, 22/52, 0, 0.05] + 
          put[s, k3, 0.2, 20/52, 0, 0.05] - 
          put[s, k2, 0.2, 20/52, 0, 0.05], {s, 40, 100}, 
        PlotRange -> All];

The Greeks (sensitivities of the option price on the parameters)

Clear[s, k, σ, T, t, r]
delta = D[call[s, k, σ, T, t, r], s] // PowerExpand // 
        Simplify

!!! Warning !!! Time consuming calculation.

theta = D[call[s, k, σ, T, t, r], t] // PowerExpand // 
        FullSimplify

rho = D[put[s, k, σ, T, t, r], r] // PowerExpand // 
        Simplify

The calculation of Vega for a CALL takes  cca 2 minutes:

vegac = D[call[s, k, σ, T, t, r], σ] // PowerExpand // 
        FullSimplify

Binomial trees

In this Section we closely follow Chapter 18, Shaw, W.: Modelling Financial Derivatives with Mathematica.

Introduction

Let us consider a time interval (t, t + Δt). Under the Black-Scholes model the conditional distribution of (Subscript[S, t + Δt]) given (Subscript[S, t]) is log-normal \(TraditionalForm\`LN(
        ln\ S\_t + \((μ - 
                1\/σ\^2)\) Δt, \(σ\^2\) Δt)\). In a risk neutral world μ = r, where r is the riskless rate and we have

E(
          (Subscript[S, t + Δt]) | (Subscript[S, t])) = (Subscript[S, t]) 
            (e^(r Δt)) = :(Subscript[E, c])     and     var(
          (Subscript[S, t + Δt]) | 
            (Subscript[S, t])) = ((Subscript[S, t])^2) (e^(2  r Δt))(
              (e^((σ^2) Δt)) - 
                1) = :(Subscript[V, c]).

Under the discrete binomial random walk we suppose that the stock price (Subscript[S, t]) moves to one of two new values, (Subscript[S, t + Δt]) = (Subscript[S, t]) u and (Subscript[S, t + Δt]) = (Subscript[S, t]) d where d < 1 < u with probabilities p and q := 1 - p, respectively. In this case we have

E(
          (Subscript[S, t + Δt]) | (Subscript[S, t])) = (Subscript[S, t])(
            p u  +  q d) = :(Subscript[E, b])    and     var(
          (Subscript[S, t + Δt]) | (Subscript[S, t])) = ((Subscript[S, t])^2) 
            p ((q(u - d))^2) = (Subscript[V, b]).

While in the continuous model we have two parameters r and (σ^2), ceteris paribus, in the discrete model we have three parameters, i.e., u,  d,  p. In practice we choose the parameters of the binomial model to match the parameters of the continuous one in terms of the moments:

(Subscript[E, b]) = (Subscript[E, c])     and     (Subscript[V, b]) = 
          (Subscript[V, c]).

The reason is simple; usually we can estimate r and (σ^2) directly which is hardly possible for the binomial model. Since there is one more degree of freedom, we have to postulate a reasonable relationship between the parameters involved. The corresponding system of equations now reads:

Clear[p, u, d, r, σ]
eqs = {p*u + (1 - p)*d == Exp[r*Δt], 
          p*(1 - p)*(((u - d))^2) == 
            Exp[2*r*Δt]*(Exp[(σ^2) Δt] - 1)};

Jarrow-Rudd scheme

Under this scheme it is supposed that both up and down movements have the same probability.

Clear[p, u, d, r, σ]
p = (1/2); solJR = Solve[eqs, {u, d}]

Only the largest solution makes sense:

{uJR, dJR} = {u, d} /. solJR[[2]]

This is an approximation:

approxJR = {Series[uJR, {Δt, 0, 1}] // Normal // 
          PowerExpand, 
Series[
              dJR, {Δt, 0, 1}] // Normal // 
          PowerExpand}

which is very close to frequently used

Exp[approxJR]/E

Cox-Ross-Rubinstein scheme

Under this scheme it is supposed that u = 1/d:

Clear[p, u, d, r, σ]
u = (1/d); 
    solCRR = Solve[eqs, {p, d}] // PowerExpand // 
        FullSimplify

After a hand-made simplification (not too difficult) we get

Clear[CRRupBin];
CRRupBin[n_, r_, q_, σ_, T_] := Module[{f = Exp[-(r-q)*T/n], 
g = Exp[(r-q+σ^2)*T/n]},
1/2*(f+g + Sqrt[(f+g)^2-4])];
Clear[CRRdownBin];
CRRdownBin[n_, r_, q_, σ_, T_] := Module[{f = Exp[-(r-q)*T/n], 
g = Exp[(r-q+σ^2)*T/n]},
1/2*(f+g - Sqrt[(f+g)^2-4])];
Discount[n_, r_, T_] := N[Exp[-r T/n]];

Binomial tree for American options

Let us consider an n steps' evolution, time to expiry T and the time step Δt = T/n. Denote (i, j) the jth node at time iΔt, i = 1,  ... , n,  j = 0,  ... , i. Let (Subscript[V, i, j]) be the value of the option (either CALL or PUT) at node (i, j). The price of the stock is (Subscript[S, 0]) (u^j) (d^(i - j)) at that node. Since the value for an American PUT is max(K - (Subscript[S, T]), 0) at expiry, the value is

(Subscript[V, n, j]) = 
        max(K - (Subscript[S, 0]) (u^j) (d^(n - j)), 0),      j = 0, 
      1,  ... , n.

With probability p the system moves from node (i, j) at time iΔt to (i + 1, j + 1) at time (i + 1) Δt and with probability 1 - p to (i + 1, j). Without possibility of an early exercise we get

(Subscript[V, i, j]) = (e^(-r Δt))(
          p (Subscript[V, i + 1, j + 1]) + (1 - p) (Subscript[V, i + 1, j])),    i = 
        0,  ... , n - 1,    j = 0,  ... , i.

If an early exercise is taken into account, the value must be compared with its intrinsic value. Thus we finally get

(Subscript[V, i, j]) = 
        max {K - (Subscript[S, 0]) (u^j) 
                (d^(i - j)), (e^(-r Δt))(
              p (Subscript[V, i + 1, j + 1]) + (1 - p) 
                  (Subscript[V, i + 1, j]))},    i = 0,  ... , 
      n - 1,    j = 0,  ... , i.

Starting from level n we go backward through the tree to the desired value (Subscript[V, 0, 0]).

Implementation

Here also a dividend yield q is considered.

$RecursionLimit=5000;

afunc[n_,  r_,  q_,  T_]  :=  
        N[Exp[(r - q) T/n]];
Pfunc[up_, down_, a_, disc_] := 
        N[disc*(a - down)/(up - down)];
Qfunc[up_, down_, a_, disc_] := 
        N[disc - Pfunc[up, down, a, disc]];

Clear[CRRAmericanOption]
CRRAmericanOption[S_, n_, σ_, T_, r_, q_, exercise_Function] :=
Module[{u = CRRupBin[n,r,q,σ,T], d = CRRdownBin[n,r,q,σ,T],
a = afunc[n,r,q,T],
disc = Discount[n,r,T],
P,Q,TreeOp,value,delta,gamma,theta},
P = Pfunc[u,d,a,disc];
Q = Qfunc[u,d,a,disc];
TreeOp[node_,level_] := TreeOp[node,level] = 
If[level==n,exercise[S d^node u^(level-node)],
Max[{P,Q}.{TreeOp[node,level+1],TreeOp[node+1,level+1]},
exercise[S d^node u^(level-node)]]];
value=TreeOp[0,0];
delta = (TreeOp[0,1] - TreeOp[1,1])/((u-d)*S);
gamma = ((TreeOp[0,2] - TreeOp[1,2])/(u^2-1) - 
(TreeOp[1,2]-TreeOp[2,2])/(1-d^2))*2/(S^2 (u^2-d^2));
theta = (TreeOp[1,2] - TreeOp[0,0])/(2*T/n);
Clear[TreeOp];
{value, delta, gamma, theta}];

CRRAmericanPut[S_, K_, n_, r_, q_, σ_, T_] := 
        CRRAmericanOption[S, n, σ, T, r, q, 
          Max[K - #, 0] &];
CRRAmericanCall[S_, K_, n_, r_, q_, σ_, T_] := 
        CRRAmericanOption[S, n, σ, T, r, q, 
          Max[# - K, 0] &];

Recall the meaning of the parameters:

CRRAmericanCall[S_,K_,n_,r_,q_,σ_,T_]
call[s_,k_,σ_,T_,t_,r_]

CRRAmericanCall[9,  10, 10,  0.1,  0,  0.2,  
      1]
call[9, 10, 0.2, 1, 0, 0.1]

Look at the convergence:

Table[{m, CRRAmericanCall[9,  10, m,  0.1,  0,  0.2,  
              1][[1]] // Evaluate}, {m, 3, 8}]

The following plot (number of steps equal to 80) shows that the convergence is not too convincing and that the recommended number of steps 30 may not be sufficient.

ListPlot[
        Table[{m, CRRAmericanCall[9,  10, m,  0.1,  0,  0.2,  
                  1][[1]] // Evaluate}, {m, 3, 80}], 
        PlotJoined -> True];

The following Example comes from HULL, p.391, where the approximate values u = (e^(σ (Sqrt[ Δ])t)), d = 1/u are used. His result (after 5 steps) is 4.49. A comparison with usual BS formula for an European CALL is also given:

CRRAmericanCall[S_,K_,n_,r_,q_,σ_,T_]
call[s_,k_,σ_,T_,t_,r_]

CRRAmericanPut[50, 50, 5, 0.1, 0, 0.4, 
      0.4167]
put[50, 50, 0.4, 0.4167, 0, 0.1]

Markowitz Mean-Variance Portfolio

We assume that N assets have expected returns given by the N×1vector r and the covariance matrix of returns V. We will seek the optimal portfolio in the sense that the Sharpe's ratio (also called Sharpe's measure of portfolio) is maximum:

max(((x^T)r)/(Sqrt[(x^T) V x]))subject to conditionsx≥0and(∑_(n = 1)^N)(Subscript[x, n])=1,

so that short sales are not allowed.

We have 8 assets traded on Prague Stock Exchange with the expected returns

returns = {-9.92488533927622*10^-04, 
          5.72524775783378*10^-03, 
          9.04818959590086*10^-03, -1.44968966764537*10^-03, -1.74045537023702*10^-02, 5.70188394543127*10^-03, 
          3.33476815091662*10^-03, 
          2.15075479228363*10^-03};

and the covariance matrix of returns

covariancematrix = {{1.35640293495777*10^-03, 
            2.56330097904098*10^-04, 3.33037161321956*10^-04, 
            1.80223328810003*10^-04, 2.81325992644196*10^-04, 
            5.93263545509886*10^-04, 7.13610711759175*10^-04, 
            3.20932535820455*10^-04}, {2.56330097904098*10^-04, 
            2.35439634554186*10^-03, 4.42407197078591*10^-04, 
            6.57567906193998*10^-06, 4.01068448708833*10^-04, 
            5.93785356343409*10^-04, 1.30961261626908*10^-04, 
            2.93338399396577*10^-04}, {3.33037161321956*10^-04, 
            4.42407197078591*10^-04, 
            1.57895651834543*10^-03, -1.16804549527263*10^-04, 
            1.85336104268171*10^-04,  .000176059157865, 
            3.41400743910472*10^-04, 
            1.86624552570071*10^-04}, {1.80223328810003*10^-04, 
            6.57567906193998*10^-06, -1.16804549527263*10^-04, 
            6.16786555694974*10^-03, 4.84828341822248*10^-04, 
            8.50333422445824*10^-04, 7.74803938060802*10^-05, 
            4.39412705238888*10^-04}, {2.81325992644196*10^-04, 
            4.01068448708833*10^-04, 1.85336104268171*10^-04, 
            4.84828341822248*10^-04, 9.45147795376653*10^-03, 
            1.78331918975426*10^-03, 5.0987274991634*10^-05, 
            1.49823427374484*10^-03}, {5.93263545509886*10^-04, 
            5.93785356343409*10^-04,  .000176059157865, 
            8.50333422445824*10^-04, 1.78331918975426*10^-03, 
            7.38754498647722*10^-03, 1.03084013803777*10^-03, 
            9.34434830622224*10^-04}, {7.13610711759175*10^-04, 
            1.30961261626908*10^-04, 3.41400743910472*10^-04, 
            7.74803938060802*10^-05, 5.0987274991634*10^-05, 
            1.03084013803777*10^-03, 2.1036055256453*10^-03, 
            2.40818723671803*10^-04}, {3.20932535820455*10^-04, 
            2.93338399396577*10^-04, 1.86624552570071*10^-04, 
            4.39412705238888*10^-04, 1.49823427374484*10^-03, 
            9.34434830622224*10^-04, 2.40818723671803*10^-04, 
            2.77726103166242*10^-03}};

We see that the above matrix is badly conditioned:

Det[covariancematrix]

Next we make use of the procedure NMaximize implemented in the Standard Numerical Package Mathematica, version 4.2. First we create the vector of unknown proportions in the portfolio:

xx = Table[(Subscript[x, i]), {i, 1, Length[returns]}]

The criterion is given by

criterion = (((xx . returns))/Sqrt)[xx . covariancematrix . xx];

This loads the package:

Needs["NumericalMath`NMinimize`"]

The constraints are

constraints1 = 
        Apply[And, 
          Table[(Subscript[x, i]) >= 0, {i, 1, 
              Length[returns]}]];
constraints2 = 
        Sum[(Subscript[x, i]), {i, 1, Length[returns]}] == 
          1;
constraints = constraints1 && (constraints2)

Next we solve the problem using the option for the method:

soldif = 
      NMaximize[{criterion, constraints}, xx, 
        Method -> "DifferentialEvolution"]

We see some problems. Nevertheless, we have got a solution satisfying the constraints:

Apply[Plus, xx /. soldif[[2]]]
Print["Expected return of the optimal portfolio = ", (xx /. 
            soldif[[2]]) . returns]

Looking on the criterion, we see that the solution is independent of a nonzero multiple. So we can suspend the second constraint and make the

constraints3 = 
      Apply[And, Table[(Subscript[x, i]) >= 0, {i, 1, Length[returns]}]]

soldif = NMaximize[{criterion, constraints3}, xx ]

We see that the optimal value is slightly higher than in the above case. The optimal portfolio together with its expected return (substantially higher than in the previous case) is simply obtained:

sol = xx /. soldif[[2]];
sol = sol/Apply[Plus, sol]
sol . returns