# Optimization of a portfolio of stocks

I would like to use Mathematica to perform an optimization of a portfolio composed by 10 stocks. I did the first part to compute yields and expected yields, but I don't know how to finish the optimization.

portf = {"AAPL", "BA", "IBM", "BMW.DE", "DIS", "R", "PEP", "BRBY.L", "AXP", "BTI"}

priceShift = FinancialData[#, "Price", {{2004}, {2011}, "Month"}][[All, 2]] & /@ portf


I am interested in the minimization of the portfolio variance (i.e. risk) for a given return.

Well, let me try to answer the OP's question. And thanks MMA.SE, for reopening this interesting question!

# DATA

To answer this question, you have to get the data using Mathematica's FinancialData function. This was the only thing originally done by the OP!

First step: define which stocks will be included in the portfolio:

Portfolio = {"AAPL", "BA", "IBM", "BMW.DE", "DIS", "R", "PEP", "BRBY.L", "AXP", "BTI"};


Second step: get the data

data = FinancialData[#, "Price", {{2004}, {2011}, "Month"}][[All, 2]] & /@ Portfolio;


You can graphically verify the price behaviour of all porfolio components for the specified period:

ListLogPlot[data, PlotRange -> All, Joined -> True] Third step: in order to perform the optimization, we have to work with returns, not prices!

Returns = Differences[Log[data[[#]]]] & /@ Range;


Note that the returns are calculated as $$Log(data[n+1])-Log(data[n])$$. This is because $$Log\left(1-\frac{data[n+1]-data[n]}{data[[n+1]]}\right)==Log\left(data[n+1]\right)-Log\left(data[n]\right)$$ and $$\frac{data[n+1]-data[n]}{data[n+1]}$$ is the arithmetic return. That is to say, what we're tracking now is not the arithmetic return, $$r_i$$, but the logarithmic return, $$Log\left(1+r_i\right)$$, which is nearly the same as $$r_i$$ for values less than one.

You can also graphically verify the returns behaviour for all portfolio components:

ListLinePlot[Returns, PlotRange -> All, ImageSize -> Medium] A second possibility would be to use BoxPlots:

BoxWhiskerChart[Returns, "Outliers", ImageSize -> Large, ChartLabels -> Portfolio] Fourth step (optional): in order to optimize (i.e., minimize) the portfolio risk, we have to find the returns relationship between all pairs of stocks. This is easily done with a correlation matrix. However, bear in mind that portfolio theory works with covariance matrices, and we will calculate them latter!

Correlation[Transpose@Returns] // MatrixForm Another possibility would be to graphically plot the correlation matrix:

Legended[MatrixPlot[Correlation[Transpose[Returns]], ColorFunction -> #], BarLegend[{#, {-1, 1}}]] &@"ThermometerColors" # THE THEORY

To perform the optimization, we have to define all necessary variables first. As we have 10 variables, we have to make the process of variables creation as simple as possible. Let's do it!

First step: define the amount of variables involved

n = 10;


Second step: we have to define a weighting vector

WeightsVector = Subscript[w, #] & /@ Range[n] Third step: we have to define a mean vector

MeanVector = Subscript[\[Mu], #] & /@ Range[n] Fourth step: we have to define a variance vector

VarianceVector = Subsuperscript[\[Sigma], #, 2] & /@ Range[n] Fifth step: we have to define a standard deviation vector

SDVector = Subscript[\[Sigma], #] & /@ Range[n] Sixth step: define the covariance matrix

CovMatrix = Array[\[Sigma], {n, n}] /. {\[Sigma][i_, j_] :> Subscript[\[Sigma], ToString@j <> ToString@i] /; i > j, \[Sigma][i_, j_] :> Subscript[\[Sigma], ToString@i <> ToString@j] /; i < j, \[Sigma][i_, j_] :>Subsuperscript[\[Sigma], i, 2] /; i == j};
CovMatrix // MatrixForm With Mathematica it's easy to verify the Markowitz equation for the portfolio risk:

PortfolioVariance := WeightsVector.CovMatrix.WeightsVector;
PortfolioVariance // Expand And the portfolio mean return can also be easily computed with Mathematica:

PortfolioMean := WeightsVector.MeanVector;
PortfolioMean # ASSIGNING VALUES TO THE VARIABLES

Step 1: assign values to the mean vector

MeanVector = Mean@Transpose@Returns*12


Please observe that multiplying the average monthly return by 12 gives the expected annual return for the stock.

Step 2: assign values to the variance vector

VarianceVector = Variance@Transpose@Returns*12


Step 3: assign values to the standard deviation vector

SDVector = StandardDeviation@Transpose@Returns*Sqrt@12


Please observe that to annualize the standard deviation you have to multiply monthly standard deviation by the square root of 12.

Step 4: assign values to the covariance matrix

CovMatrix = Covariance@Transpose@Returns


# OPTIMIZATION

Before doing the optimization, we can use Mathematica to verify the Markowitz equation, now with all assigned values for the variables:

PortfolioVariance := WeightsVector.CovMatrix.WeightsVector;
PortfolioVariance // Expand As we see, the only thing we have to find is exactly the weight that will be assigned for every portfolio component. We know that the weights of all portfolio components must sum up to one, and we will use it as a constraint for our optimization problem:

MV = FindMinimum[{PortfolioVariance, Total[WeightsVector] == 1}, WeightsVector] Finally, in order to better interprete the results, we can plot the weights using a BarChart:

BarChart[{Last[MV[[2, #]]] & /@ Range[n]}, ImageSize -> Large, ChartLabels -> Portfolio] As we see, in order to optimize the portfolio risk we have to buy 6 out of the 10 stocks and short sell 4 of them (those stocks with negative weights).

• Rod, this is a really nice answer, I am glad that this question is reopened: this thread could be the flagship of further finance-related questions. Two minor advices: your variables should start with lowercase letters (Returns is really close to Return). Second: try to be more concise. No need to denote each code line with a comment or with a title. The more concise and fluid your code is, the easier it is to copy and paste it by others. And n=10 really doesn't need comments. Dec 12, 2013 at 18:23
• @IstvánZachar Thanks István! My intention was actually to write my answer in an "paper"-like format... that's why I explained every line of code. But I do agree with you: sometimes I need to be more "concise"... :-)
– Rod
Dec 12, 2013 at 20:49
• this answer is outstanding May 11, 2015 at 22:44
• Changing "2011" to "2016" seems to error now... any idea why? Jan 1, 2017 at 0:28
• @ István Zachar - As someone who's trying to follow the logic, I strongly disagree with your suggestion about being concise. This is heavy stuff and I would rather have someone like Rod overestimate the need for comments than try to muddle through a big block of minimally commented code. Sep 17, 2020 at 19:34
portf = {"AAPL", "BA", "IBM", "BMW.DE", "DIS", "R", "PEP", "BRBY.L", "AXP", "BTI"};
prices = FinancialData[#, "Price", {{2004}, {2011}, "Month"}] & /@  portf;


Returns are often calculated as the difference of the logarithms of the prices:

Differences[Log@#] &[prices[[1, All, 2]]]


This works because $\log(\text{price}_{new})-\log(\text{price}_{old})=\log(\text{price}_{old}(1+r))-\log(\text{price}_{old})=\log(1+r)$, which is

Series[Log[o (1 + r)] - Log[o], {r, 0, 2}] that is, approximately $r$.

But we can get returns from FiancialData right away using "FractionalChange" instead of "Price":

portReturns = FinancialData[#, "FractionalChange", {{2004}, {2011}, "Month"}] & /@ portf;


Can't hurt to check the data a bit. Let's see if we have data from the same dates:

Outer[
Length[Intersection[#1, #2]] &,
portReturns[[All, All, 1]],
portReturns[[All, All, 1]],
1
] // MatrixForm It appears not all periods are exactly the same. So, if we want to be prudent we could use only the data for the periods that match.

overlappingDates = Intersection @@ (#[[All, 1]] & /@ portReturns);
portReturnsCleaned = Select[#, MemberQ[overlappingDates, #[]] &][[All, 2]] & /@ portReturns;


Now, introduce symbolic weights and define portfolio return and variance in terms of them. The portfolio variance is the sum of the covariance matrix elements.

weights = {w1, w2, w3, w4, w5, w6, w7, w8, w9, w10};
portAverageReturns = 12 weights (Mean /@ portReturnsCleaned) // Total
portVar = 12 Total[Outer[Times, weights, weights] Covariance[portReturnsCleaned\[Transpose]], 2]


Let's find the lowest variance for a bunch of given returns (the so-called efficient frontier). The range of this will lie between:

12 Min@(Mean /@ portReturnsCleaned)


0.08960473674

12 Max@(Mean /@ portReturnsCleaned)


0.7151653802

(should have bought Apple stock)

The no short selling, no borrowing/lending case (weights sum to 1 and each weight is between 0 and 1):

effFrontier =
Table[
{Sqrt[
NMinimize[
{portVar,
Total[weights] == 1 &&
portAverageReturns == i &&
And @@ (0 <= # <= 1 & /@ weights)
}, weights
][]
], i},
{i, 0.09, 0.7, .02}] // Quiet


I collected the square root of the variance (the standard deviations) here as this is usually considered as the measure of risk.

Find the minimum variance portfolio:

minimumVariance1 = {Sqrt[portVar], portAverageReturns} /.
NMinimize[{portVar, Total[weights] == 1 && And @@ (0 <= # <= 1 & /@ weights)}, weights][]


{0.1292044085, 0.127437124}

Plotting this with 10,000 randomly selected portfolios:

Show[
ListLinePlot[effFrontier,
Epilog -> {Red, PointSize -> Large, Point[minimumVariance]},
Frame -> True,
FrameLabel->{"σ","Expected portfolio return","No short selling, no borrowing/lending",""},
PlotRange -> {{0, Automatic}, Automatic}, BaseStyle -> 16,
ImageSize -> 500],
ListPlot[
Table[{Sqrt[portVar], portAverageReturns} /.
Thread[weights -> #/Total[#] &[RandomReal[{0, 1}, 10]]], {10000}]]
] It appears pretty difficult to blindly guess an optimal portfolio.

Dropping the constraint of no short selling (positive weights) gets us better results:

effFrontier2 =
Table[{Sqrt[
NMinimize[{portVar, Total[weights] == 1 && portAverageReturns == i}, weights][[
1]]], i},
{i, 0.09, 0.7, .02}] // Quiet;
minimumVariance2 = {Sqrt[portVar], portAverageReturns} /.
NMinimize[{portVar, Total[weights] == 1}, weights][] One more step is to introduce lending/borrowing at the risk free rate:

portAverageReturns2 = α rf + 12 Total[(1 - α) weights (Mean /@ portReturnsCleaned)]
portVar2 = 12 Total[Outer[Times, (1 - α) weights, (1 - α) weights]
Covariance[portReturnsCleaned\[Transpose]], 2]


Assuming a risk free rate of 1%:

Block[{rf = .01},
effFrontier3 =
Table[{Sqrt[
NMinimize[{portVar2,
Total[weights] == 1 && portAverageReturns2 == i &&
And @@ (0 <= # <= 1 & /@ weights)},
Join[weights, {α}]][]], i}, {i, 0.01, 0.7, .05}];
effFrontier4 =
Table[{Sqrt[
NMinimize[{portVar2,
Total[weights] == 1 && portAverageReturns2 == i &&
And @@ (# <= 1 & /@ weights)}, Join[weights, {α}]][[
1]]], i}, {i, 0.01, 0.7, .05}];
]

Show[
ListLinePlot[effFrontier,
Epilog -> {Red, PointSize -> Large, Point[minimumVariance],
Point[minimumVariance2]}, Frame -> True, PlotStyle -> Dashed,
PlotRange -> {{0, Automatic}, Automatic}, BaseStyle -> 16,
ImageSize -> 500,
FrameLabel -> {"σ", "Expected portfolio return", "", ""},
PlotLegends ->
Placed[LineLegend[{Blue, Directive[Blue, Dashed], Darker@Green,
Directive[Darker@Green,
Dashed]}, {"Short & no borrowing/lending",
"No short & no borrowing/lending",
"Short & with borrowing/lending",
"No short & with borrowing/lending"}, LabelStyle -> 12],
Scaled[{0.25, 0.8}]]],
ListLinePlot[effFrontier2, Frame -> True,
PlotRange -> {{0, Automatic}, Automatic}, BaseStyle -> 16,
ImageSize -> 500],
ListLinePlot[effFrontier3, Frame -> True,
PlotRange -> {{0, Automatic}, Automatic}, BaseStyle -> 16,
ImageSize -> 500, PlotStyle -> Directive[Green, Dashed]],
ListLinePlot[effFrontier4, Frame -> True,
PlotRange -> {{0, Automatic}, Automatic}, BaseStyle -> 16,
ImageSize -> 500, PlotStyle -> Green],
ListPlot[Table[{Sqrt[portVar], portAverageReturns} /.
Thread[weights -> #/Total[#] &[RandomReal[{0, 1}, 10]]], {10000}]]
] • I was on vacations and couldn't see your answer before (nice one, btw). However, I'm still wondering why short sales have changed your efficient frontier if borrowing/lending is not allowed. I might be wrong, but as far as I remember there is no difference in the efficient frontier without borrowing/lending.
– Rod
Jan 7, 2014 at 11:56
• @rod That may be the case for the two-asset portfolio, which is what is generally plotted in textbooks. Generally, this is not true. Check these three sources: 1 (first picture), 2 (section 3.iii), 3 (page 10). Jan 19, 2014 at 0:36
• A note on effFrontier3 and 4, with a risk-free rate. weights limits the risky allocations but not a, the risk-free asset. It can turn >1 for portf mean returns <0.01, e.g. 0.005 which means shorting the tangency portfolio and buying the RF. Setting it to Join[weights, {α}] in the restrictions seems to work, apart from some accuracy error messages. I know the mathematical solution but preventing short-selling in the stocks but not in the whole portfolio still looks counter-intuitive. Am I getting something wrong here? And congratulations for an extremely useful answer. May 31, 2019 at 1:38
• It appears not all periods are exactly the same. So, if we want to be prudent we could use only the data for the periods that match. The reason that the periods aren't exactly the same is because two of the issues, BMW.DE and BRBY.L, are foreign stocks and have a different holiday calendar. There is no justification for removing the returns for dates that don't align. We are interested in minimizing the variance of the stocks and you're deleting important information here. If we had a huge swing of AAPL during a British Holiday, your portfolio wouldn't reflect that added risk. Sep 24, 2020 at 22:24

This is nothing more than a fusion of the excellent answers by Rod and Sjoerd C. de Vries updated for MMa 12.

Extract the prices from Wolfram's data service.

portfolio = {"AAPL", "BA", "IBM", "BMW.DE", "DIS", "R", "PEP", "BRBY.L", "AXP", "BTI"};
assetCount = Length[portfolio];
priceHistory = (QuantityMagnitude[FinancialData[#1, "Price", {{2016}, {2020}, "Month"}]["Path"][[All,2]]] & ) /@ portfolio;
ListLogPlot[priceHistory, PlotRange -> All, Joined -> True] Calculate the actual returns based on the log of the difference in the asset price from month to month.

actualReturns = (Differences[Log[priceHistory[[#1]]]] & ) /@ Range[assetCount];
ListLinePlot[actualReturns, PlotRange -> All] Note that what we're tracking now is not the arithmetic returns, but the logarithmic returns. $$actualReturns = Log_{10}(1+r_i)$$ Which is close to $$r_i$$ when $$r_i < 1$$ and approximately $$Log_{10}(r_i)$$ when the value is large, so it acts as a normalization function for the returns.

Also note that I disagree strongly with Sjoerd C. de Vries here. He's deleting valuable information about the variance of the assets with his attempt to 'clean up' the data. If AAPL has a wild swing on during a British holiday, his cleaning step will make AAPL appear less volatile than it really is.

The variables used for the asset's weights.

weights = Table[Subscript[w, i], {i, Range[assetCount]}, {j, 1}];


This is the formula used to minimize the mean variance.

covariance = Covariance[Transpose[actualReturns]];
portfolioVariance = Transpose[weights] . covariance . weights
portfolioVariance = Expand[portfolioVariance[][]];


Note that I've set up the vectors as column vectors which is what most textbooks do when introducing this subject. What we're actually doing here is this: $$\mathbf w'\mathbf\Sigma \mathbf w$$

Where $$\mathbf w$$ is the vector of weights and $$\mathbf\Sigma$$ is the covariance matrix.

Minimum and Maximum Returns

In order to chart the returns, we need to know the range of possible returns. The minimum and maximum return for a given portfolio are simply the best and worst annual return of an asset while setting the weight of all the other assets to zero.

minimumReturn = 12*Min[Mean /@ actualReturns];
maximumReturn = 12*Max[Mean /@ actualReturns];


Calculate the Efficient Frontier.

yearlyWeightedReturns = First[Total[12*weights*Mean /@ actualReturns]];
efficientFrontier = Table[{Sqrt[FindMinimum[{portfolioVariance, First[Total[weights]] == 1 && yearlyWeightedReturns == i}, weights][]], i}, {i, minimumReturn, maximumReturn, 0.02}];


Now let's see how changing the desired return effects the risk.

Manipulate[minimum = FindMinimum[{portfolioVariance, Total[weights] == 1 && yearlyWeightedReturns == return}, weights];
plot1 = BarChart[{(Last[minimum][[#1,2]] & ) /@ Range[assetCount]}, PlotRange -> {-0.4, 0.6}, AxesOrigin -> {0, -0.4}, ImageSize -> Large, ChartLabels -> portfolio];
minimumVariance = First[minimum]; plot2 = ListLinePlot[efficientFrontier, Epilog -> {Red, PointSize -> Large, Point[{Sqrt[First[minimum]], return}]}, Frame -> True,
FrameLabel -> {"\[Sigma]", "Expected portfolio return"}, PlotRange -> {{0, Automatic}, Automatic}, BaseStyle -> 16, ImageSize -> Large]; GraphicsGrid[{{plot1}, {plot2}}],
{{return, 0.1, "Desired Return"}, minimumReturn, maximumReturn, Appearance -> "Labeled"}] 