## Overview

Portfolio Theory was introduced by Harry Markowitz in 1952, which later awarded him a Nobel Prize in economics. Portfolio theory is a neat application of Lagrange multipliers. The objective function is second order and all constraints are first order. After applying Lagrange multiplier, it becomes a system of linear equations.

## Problem Statement

Given N random variables with known expected value u[i],i=1..N and covariance matrix v[i,j] for i,j=1..N and a target return R. Find weight w[i] such that

1. sum[i=1..N]( w[i] ) = W

2. sum[i=1..N]( w[i]*u[i] ) = R

3. minimize sum[i=1..N,j=1..N]( w[i]*v[i,j]*w[j] ) i.e. variance of portfolio

(W is sum of weight, usually set W = 1)

## Solution

let L(w,p1,p2) = sum[i=1..N,j=1..N]( w[i]*v[i,j]*w[j] )

+ p1*(sum[i=1..N]( w[i] ) – W)

+ p2*(sum[i=1..N]( w[i]*u[i] ) – R)

Differentiate L(w,p1,p2) w.r.t w[i] and set it to zero, we have

2*v[i,i]*w[i] + sum[j=1..N,i!=j]( w[j]*v[j,i] + v[i,j]*w[j] ) + p1 + p2*u[i] = 0

Since covariance matrix is symmetric, simplify, we have

sum[j=1..N]( 2*v[i,j]*w[j] ) + p1 + p2*u[i] = 0 for i=1..N — equation1

together with the 2 original constraints

1. sum[i=1..N]( w[i] ) = W

2. sum[i=1..N]( w[i]*u[i] ) = R

We have a N+2 independent linear equations for N+2 variables, so we have a system of linear function.

To find the assignment of w[i], we form the matrix of size (N+2)*(N+2) s.t.

M[i,j] = 2*v[i,j] for i,j=1..N

M[N+1,i] = M[i,N+1] = 1 for i=1..N

M[N+2,i] = M[i,N+2] = u[i] for i=1..N

else 0

(NOTE: this matrix is positive semidefinite)

And then a vector vec s.t.

vec[N+1] = W

vec[N+2] = R

else 0

And then solve inverse(M)*vec for answer. Use any decent linear solver for the job, usually solve in O(N^2)~O(N^3).

## Code

import numpy as np
def portfolio(u,v,R,W=1):
n = len(u)
m = np.zeros((n+2,n+2))
m[0:n,0:n] = 2*v
m[n,0:n] = m[0:n,n] = 1
m[n+1,0:n] = m[0:n,n+1] = u
v = np.zeros(n+2)
v[n] = W
v[n+1] = R
return np.linalg.solve(m, v)

## Special Case

If all random variables are independent, then v[i,j] = 0 for i != j, we could find the solution in O(N) without matrix inverse.

Denote v[i] = v[i,i]

Simplify equation1, we have

2*v[i]*w[i] + p1 + p2*u[i] = 0

=> w[i] = -(p1+p2*u[i])/(2*v[i]) — equation2

(if some v[i] == 0, this method does not work)

put this into 1 & 2, we have

1. sum[i=1..N]( -(p1+p2*u[i])/(2*v[i]) ) = W

2. sum[i=1..N]( -(p1*u[i]+p2*u[i]^2)/(2*v[i]) ) = R

simplify, we have

let A = sum[i=1..N]( 1/v[i] )

let B = sum[i=1..N]( u[i]/v[i] )

let C = sum[i=1..N]( u[i]^2/v[i] )

1. A*p1 + B*p2 = -2*W

2. B*p1 + C*p2 = -2*R

let D = A*C – B^2, then

p1 = -2*(C*W-B*R)/D

p2 = -2*(A*R-B*W)/D

put them back to equation2

let t1 = C*W-B*R

let t2 = A*R-B*W

w[i] = (t1+t2*u[i])/D/v[i]

### Code

def portfolioIndependent(u,v,R,W=1):
inv = 1./v
a = inv.sum()
b = (u*inv).sum()
c = (u*u*inv).sum()
d = a*c - b*b
t1 = c*W - b*R
t2 = a*R - b*W
return (t1+t2*u)/d/v