Quadratic forms, namely polynomials whose terms are all of degree of two, are of great importance in linear algebra and its various applications. In matrix notation, a quadratic form for a \(K\)-dimension vector \(\vec{w} = \t{\left( w_1, w_2, \dots, w_K \right)}\) can be written as:
\[ \fun{q}{\vec{w}} = \t{\vec{w}} \vec{A} \vec{w} \] where the matrix \(\vec{A}\) specifies the coefficients in the quadratic form. When expanded, we can clearly see how the above expression corresponds to a quadratic form:
\[ \begin{align*} \fun{q}{\vec{w}} &= \t{\vec{w}} \vec{A} \vec{w} \\ &= \t{\begin{pmatrix} w_1 & w_2 & \dots & w_K \end{pmatrix}} \begin{pmatrix} a_{11} & a_{12} & \cdots & a_{1K} \\ a_{21} & a_{22} & \cdots & a_{2K} \\ \vdots & \vdots & \vdots & \vdots \\ a_{K1} & a_{K2} & \cdots & a_{KK} \\ \end{pmatrix} \begin{pmatrix} w_1 \\ w_2 \\ \vdots \\ w_K \\ \end{pmatrix} \\ &= \t{\begin{pmatrix} w_1 & w_2 & \dots & w_K \end{pmatrix}} \begin{pmatrix} a_{11} w_1 + a_{12} w_2 + \cdots + a_{1K} w_K \\ a_{21} w_1 + a_{22} w_2 + \cdots + a_{2K} w_K \\ \vdots \\ a_{K1} w_1 + a_{K2} w_2 + \cdots + a_{KK} w_K \\ \end{pmatrix} \\ &= \begin{matrix} & a_{11} w_1 w_1 &+& a_{12} w_1 w_2 &+& \cdots &+& a_{1K} w_1 w_K \\ + & a_{21} w_2 w_1 &+& a_{22} w_2 w_2 &+& \cdots &+& a_{2K} w_2 w_K \\ + & \dots &+& \dots &+& \dots &+& \dots \\ + & a_{K1} w_K w_1 &+& a_{K2} w_K w_2 &+& \cdots &+& a_{KK} w_K w_K \end{matrix} \end{align*} \]
A quadratic form becomes particular interesting when its coefficient matrix \(\vec{A}\) is a covariance matrix — suppose we have a random vector \(\vec{X} = \t{\left( X_1, X_2, \dots, X_K \right)}\), and let \(\vec{Y}\) be a linear combination of the elements in \(\vec{X}\) weighted by \(\vec{w}\), then \(\vec{Y}\) is also a random variable with the following expectation and variance: \[ \begin{align*} \expect{\vec{Y}} &= \expect{\t{\vec{w}} \vec{X}} \\ &= \t{\vec{w}} \expect{\vec{X}} \\ \vari{\vec{X}} &= w_1 w_1 \covari{X_1, X_1} + w_1 w_2 \covari{X_1, X_2} + \dots + w_K w_K \covari{X_K, X_K} \\ &= \begin{pmatrix} w_1 \covari{X_1, X_1} & w_1 \covari{X_1, X_2} & \cdots & w_1 \covari{X_1, X_K} \\ + & + & \ddots & + \\ w_2 \covari{X_2, X_1} & w_2 \covari{X_2, X_2} & \cdots & w_2 \covari{X_2, X_K} \\ + & + & \ddots & + \\ w_K \covari{X_K, X_1} & w_K \covari{X_K, X_2} & \cdots & w_K \covari{X_K, X_K} \end{pmatrix} \begin{pmatrix} w_1 \\ w_2 \\ \vdots \\ w_K \end{pmatrix} \\ &= \begin{pmatrix} w_1 & w_2 & \cdots & w_K \end{pmatrix} \begin{pmatrix} \covari{X_1, X_1} & \covari{X_1, X_2} & \cdots & \covari{X_1, X_K} \\ \covari{X_2, X_1} & \covari{X_2, X_2} & \cdots & \covari{X_2, X_K} \\ \vdots & \vdots & \ddots & \vdots \\ \covari{X_K, X_1} & \covari{X_K, X_2} & \cdots & \covari{X_K, X_K} \end{pmatrix} \begin{pmatrix} w_1 \\ w_2 \\ \vdots \\ w_K \end{pmatrix} \\ &= \t{\vec{w}} \vec{\Sigma} \vec{w} \end{align*} \]
In other words, the quadratic form \(\t{\vec{w}} \vec{\Sigma} \vec{w}\) represents the variance of a linear combination of a random vector’s constituents when its covariance matrix is \(\vec{\Sigma}\).
This representation of variance of a linear combination of random variables often comes up in the context of financial portfolio optimisation. The usual set up is that we have a collection of \(K\) assets whose returns over a unit period of time, as defined by \(R_t = \frac{P_{t} - P_{t-1}}{P_{t-1}}\), are represented by the elements in a random vector \(\vec{X} = \left(X_1, X_2, \dots, X_K\right)\). A portfolio can be formed from these assets by allocating a weight to each aseet according to a \(K\)-dimensional weight vector \(\vec{w}\); the return of the portfolio is therefore a linear combination of the elements of \(\vec{X}\), i.e. \(\t{\vec{w}} \vec{X}\); the variance of the portfolio’s return, as we have shown above, is exactly \(\t{\vec{w}} \vec{\Sigma} \vec{w}\).
Once the expectation and the variance of the portfolio are defined, we can do some interesting things with them. The variance/standard-deviation of the return, for example, is often used as a measure for the risk associated with the portfolio — a desriable portfolio would have a high expected return with a low level exposure to risk. In its simplest form, this can be translated into the following optimisation problem:
\[ \begin{array}{rll} \text{maximise: } & \qquad \frac{1}{2} \t{\vec{w}} \vec{\Sigma} \vec{w} & & \text{(we want to minimise the variance of the portfolio)} \\ \text{subject to: } & \qquad \t{\vec{1}} \vec{w} &= 1 & \text{(this is to ensure the weights sum up to a constant)} \\ & \qquad \t{\expect{\vec{R}}} \vec{w} &\ge C & \text{(the expected return of the portfolio should be greater than a threshold)} \end{array} \]
Note that here we are not considering any risk-free assets and the transaction costs incurred when trading into the desired position. This particular formulation of the portfolio optimisation can be solved with quadratic programming. To demonstrate the process, I will replicate part of the example from Matlab’s documentation on portfolio optimisation with R language.
I trasncribed the original dataset used in the Matlab documentation into a RData file and load it here so that we can reproduce some of the results:
load('BlueChipStockMoments.RData')Just take a quick look at the variables loaded from the RData file. The dataset contains the means and the covariance matrix of the returns of the 30 assets:
print(asset_list) # print the list of tickers for the assets## [1] "AA" "AIG" "AXP" "BA" "C" "CAT" "DD" "DIS" "GE" "GM"
## [11] "HD" "HON" "HPQ" "IBM" "INTC" "JNJ" "JPM" "KO" "MCD" "MMM"
## [21] "MO" "MRK" "MSFT" "PFE" "PG" "SBC" "UTX" "VZ" "WMT" "XOM"print(asset_means)## AA AIG AXP BA C
## 0.0149243571 0.0095278690 0.0113101190 0.0069712619 0.0137759762
## CAT DD DIS GE GM
## 0.0150674167 0.0029306667 0.0030687738 0.0090499048 0.0058126786
## HD HON HPQ IBM INTC
## 0.0144927857 0.0079065714 0.0079954286 0.0130997143 0.0146523810
## JNJ JPM KO MCD MMM
## 0.0112448452 0.0084923333 -0.0012911071 0.0078821786 0.0123067738
## MO MRK MSFT PFE PG
## 0.0129688214 0.0007541071 0.0150997738 0.0045406667 0.0082302738
## SBC UTX VZ WMT XOM
## 0.0033750476 0.0179275833 0.0054528095 0.0156251786 0.0093157619print(asset_covar[1:5, 1:5])## AA AIG AXP BA C
## AA 0.013675787 0.002489241 0.003981606 0.004091740 0.004833671
## AIG 0.002489241 0.005999770 0.004255254 0.002388906 0.004387941
## AXP 0.003981606 0.004255254 0.006612873 0.004149297 0.005536310
## BA 0.004091740 0.002388906 0.004149297 0.009567257 0.002800733
## C 0.004833671 0.004387941 0.005536310 0.002800733 0.008938242Two additional variables, cash_mean and cash_variance, capture the expected return and the risk associated with holding cash, a relatively riskless asset:
print(cash_mean)## [1] 0.002734881print(cash_variance)## [1] 2.494875e-06We will first approach the portfolio optimisation problem by minimising the variance of the portfolio while keeping the expected return of the portfolio above the risk-free asset (i.e. cash). We will be using the solve.QP function in the quadprog library to solve this:
pacman::p_load("tidyverse")
pacman::p_load("quadprog") # we will be using the `quadprog` library to solve the optimisation problemThe documentation of the solve.QP function is shown below:
- The
Dmatargument is simply the covariance matrix of the assets,asset_covar. - We will not be using the
dvecargument since we don’t have the first-order term in our optimisation problem, so we will fill a vector of zeros here. - The
Amatand thebvecarguments are used to specify the constraints. Some extra bits of care are needed here:- We will encode the equality constraint, \(\t{\vec{1}} \vec{w} = 1\), using the first row of
Amatand the first element ofbvec; the summation can be conveniently coded as a dot product with a vector of ones. It is also important to set themeqargument to1to indicate that the first contraint is an equality constraint. - We use the second row of
Amatto encode the inequality contraint that the expected return of the portfolio should exceed that of the risk-free asset, namelycash_mean. - The rest of the
Amatmatrix encodes the extra constraint that no shorting is allowed — without such constraint, it is possible for the solver to return negative weights for certain assets; this would mean that the portfolio will have to enter a short position on those assets by borrowing them from the market. The constraint is specified by appending an identity matrix toAmatand a vector of zeros tobvec.
- We will encode the equality constraint, \(\t{\vec{1}} \vec{w} = 1\), using the first row of
The specification of the optimisation problem then becomes quite straightforward:
K <- length(asset_means)
Amat <- cbind(1, asset_means, diag(rep(1, K)))
bvec <- c(1, cash_mean, rep(0, K) )
portfolio <- solve.QP(Dmat = asset_covar, dvec = rep(0, K), Amat = Amat, bvec = bvec, meq = 1)
print(portfolio)## $solution
## [1] 3.412747e-17 5.370080e-18 -8.939030e-17 -6.564778e-18 5.887635e-17
## [6] 4.141009e-18 -4.355909e-18 4.218599e-02 5.050503e-02 -8.400927e-19
## [11] -2.131723e-19 3.841085e-17 1.690746e-02 -2.737338e-19 7.461713e-04
## [16] 1.934536e-02 -2.077243e-17 -1.773789e-17 -2.803253e-17 8.028406e-02
## [21] 7.652513e-02 -1.583827e-19 2.344458e-02 1.166485e-01 1.911548e-01
## [26] 1.379132e-18 -1.671579e-16 0.000000e+00 6.425835e-02 3.179946e-01
##
## $value
## [1] 0.0006057603
##
## $unconstrained.solution
## [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
##
## $iterations
## [1] 20 0
##
## $Lagrangian
## [1] 1.211521e-03 0.000000e+00 7.232438e-04 2.319600e-04 7.785867e-04
## [6] 4.860429e-04 6.688958e-04 4.764734e-04 2.613163e-04 0.000000e+00
## [11] 0.000000e+00 2.538115e-04 2.706206e-05 1.080283e-03 0.000000e+00
## [16] 2.429306e-04 0.000000e+00 0.000000e+00 8.471320e-04 1.899856e-04
## [21] 2.761950e-04 0.000000e+00 0.000000e+00 3.388052e-05 0.000000e+00
## [26] 0.000000e+00 0.000000e+00 1.029870e-04 6.155201e-04 7.824339e-06
## [31] 0.000000e+00 0.000000e+00
##
## $iact
## [1] 1 5 7 29 14 3 6 19 8 9 16 21 4 20 12 28 24 13 30The portfolio$solution element captures the weights for the assets in the portfolio. In the original example from Matlab’s documentation, instead of finding a single portfolio that exhibits the lowest variance, a collection of portfolios were created, each for a different expected return level. These portfolios could then be used to plot the efficient frontier for the collection of assets so that one can see the relationship between the risks and expected returns of the various portfolios. A Matlab program equivalent to the above R program is listed as follows:
load BlueChipStockMoments.mat
p = Portfolio('AssetList', AssetList, 'RiskFreeRate', CashMean);
p = setAssetMoments(p, AssetMean, AssetCovar);
p = setDefaultConstraints(p);
% note when we ask for a single portfolio from the frontier estimation procedure,
% it will return the portfolio with minimal variance:
pwgt = estimateFrontier(p, 1);
pwgtwith the following output:
pwgt =
0
0
0
0
0
0
0
0.0422
0.0505
0
0
0
0.0169
0
0.0007
0.0193
0
0
0
0.0803
0.0765
0
0.0234
0.1166
0.1912
0
0
0
0.0643
0.3180We can see that the weights estimated using Matlab are consistent with what we obtained with our R program:
as.matrix(round(portfolio$solution, digits = 4), ncol = 1)## [,1]
## [1,] 0.0000
## [2,] 0.0000
## [3,] 0.0000
## [4,] 0.0000
## [5,] 0.0000
## [6,] 0.0000
## [7,] 0.0000
## [8,] 0.0422
## [9,] 0.0505
## [10,] 0.0000
## [11,] 0.0000
## [12,] 0.0000
## [13,] 0.0169
## [14,] 0.0000
## [15,] 0.0007
## [16,] 0.0193
## [17,] 0.0000
## [18,] 0.0000
## [19,] 0.0000
## [20,] 0.0803
## [21,] 0.0765
## [22,] 0.0000
## [23,] 0.0234
## [24,] 0.1166
## [25,] 0.1912
## [26,] 0.0000
## [27,] 0.0000
## [28,] 0.0000
## [29,] 0.0643
## [30,] 0.3180