# Kelly Criterion for Optimal Portfolio Weights with No Short Sales in Mathematica

Kelly Criterion for Optimal Portfolio Weights with No Short Sales

Summary

1. The U.S. Healthcare index seeks to track the investment results of healthcare providers, biotechnology companies, and manufacturers of advanced medical devices and pharmaceuticals.

2. GROWTH of 10,000 USD SINCE 31 JULY, 2007. Introduction

This notebook shows the relationship between Markowitz portfolio optimization and Kelly optimization. It is based on a longer whitepaper by Dr Thomas Starke (https://plot.ly/ipython-notebooks/markowitz-portfolio-optimization/) All calculation are written in Mathematica 10.

There are two sections of the demonstration; random data and actual market data. Random data; that will help you to get a idea of how to use this investment strategy for asset allocation and to improve the understanding of the modelling. Then, we shall create a backtest that the asset-weights in a kelly-optimization way.

Simulations

Setup Data

Assume that we have 4 assets, each with a return series of 500 observations. We can use NormalDistribution[] to sample returns from a normal distribution.    Random Weights and the Returns/Risk of Each Portfolio

We can generate a wide range of random weights, which can be used to create different portfolios with a wide range of returns(mean) and risks(standard deviation). We assumed that all capital shall be invested, the sum of the weights is equal to 1.    Evaluate how many of these random portfolios would perform

In this section, the goal we are calculating the mean returns as well as the standard deviation (volatility). You can also see that there is a filter that only allows to plot portfolios with a standard deviation of < 1.5 for better illustration. Generating the mean returns and volatility for 500 random portfolios:  The Efficient Frontier

S.T.    ∑w = 1
w.avgReturns = tgtReturn

where vectorOfWeights is a vector of weights, mS is the sample variance/covariance matrix, optimalReturns is a vector of returns, and muRange is the target return.

The problem is that FindMinimum gives slightly different weights than NMinimize. Both give slightly different weights; I would prefer to use FindMinimum as it is much faster, and it seems to be more accurate. Optimal Portfolio Weights   Backtesting on real market data

Retrieve data from Yahoo Finance

Optimal Portfolio Weights

Assume that we have 20 assets(U.S. Healthcare Sector), each with a return series of  3500 calendar days (2411 trade-days). We can use FinancialData to sample data from a wolfram server.   Daily Return of Assets  Kelly Criterion:

Initialization PortfolioWeight,  every 70 trade-days update and rebalance its portfolio in a Markowitz-optimal way.       Conclusions

Even through the 2008 financial crisis, The volatility and return are pretty stable.  This is most likely due to the universe selection and Increasing the number of stocks in the universe might reduce the volatility as well.

# Cholesky Decomposition

Bayesian Methods for Empirical Macroeconomics
A Course Prepared the Global Initiative Program, Bank of Korea Academy
http://personal.strath.ac.uk/gary.koop/BoK_course.html

```
#r @"../packages/MathNet.Numerics.3.2.3/lib/Net40/MathNet.Numerics.dll"
#r @"../packages/FSPowerPack.Core.Community.3.0.0.0/Lib/Net40/FSharp.PowerPack.dll"

open System
open System.Collections.Generic
open System.Data
open MathNet.Numerics.Distributions

// http://www.quantstart.com/articles/Cholesky-Decomposition-in-Python-and-NumPy
// http://rosettacode.org/wiki/Cholesky_decomposition#OCaml
// http://www.math.sjsu.edu/~foster/m143m/cholesky.pdf
// http://stackoverflow.com/questions/8565489/improvement-of-matrix-calculation-in-f

// Algorithm
let cholesky (matA:matrix) =
let x = matA.NumRows
let mutable s = 0.0
let res = Matrix.create x x 0.0
for i in 0 .. x - 1 do
for j in 0 .. i do
s <- 0.0
for k in 0 .. j - 1 do
s <- s + (res.[i,k] * res.[j,k])
if i = j then res.[i,j] <- sqrt(matA.[i,i] - s)
else res.[i,j] <- ((matA.[i, j] - s) / res.[j,j])
res

//let a =  matrix [  [25.0; 15.0; -5.0];
//                    [15.0; 18.0;  0.0];
//                    [-5.0;  0.0; 11.0] ];

let a = matrix [    [16.0; 4.0;  4.0;  -4.0];
[4.0; 10.0;  4.0;  2.0];
[4.0; 4.0; 6.0; -2.0];
[-4.0; 2.0; -2.0; 4.0] ];

let abc = cholesky a;;

```

# Numerical Integration (Monte Carlo vs Trapezoidal’s Rule)

This is the exercise 2 of question 1, the details please refer the web site of Professor Gary Koop. I found that I have done something wrong, the answer of Monte Carlo is different with Trapezoidal’s Rule. If anyone else find the solution, please feel free to leave me a message.

Bayesian Methods for Empirical Macroeconomics
A Course Prepared the Global Initiative Program, Bank of Korea Academy
http://personal.strath.ac.uk/gary.koop/BoK_course.html

```(* ---------------------- *)
(* Question 1 Ex2  *)
(* -------------------------- *)
#r "System.Windows.Forms.DataVisualization.dll"
#r @"../packages/MathNet.Numerics.3.2.1/lib/Net40/MathNet.Numerics.dll"

open System
open System.Net
open System.Windows.Forms
open System.Windows.Forms.DataVisualization.Charting
open Microsoft.FSharp.Control.WebExtensions
open MathNet.Numerics.Distributions
open MathNet.Numerics.Statistics

let theta_mean = 1.0
let theta_var = 4.0

(* the number of draws we obtain will affect the precision *)
(* of the approximation using Monte Carlo integration *)
let number_of_draws = 1000000

//f : function
//a : lower limit
//b : upper limit
//n : count of intervals
let Integrate (f : float -> float)  (a:float) (b:float) (n:int) =
let h = (b - a) / float n //step width
let TrapezoidalRule (f: float -> float) (tp:float*float) =
((snd tp) - (fst tp)) * ((f (fst tp) + f (snd tp)) / 2.0)
let startPoints = [a..h..b-h]
let endPoints = [a+h..h..b]
let intervals = List.zip startPoints endPoints
List.map (TrapezoidalRule f) intervals |> List.sum

// (x.^2).*(0.1995*exp((-(x-1).^2)/8))
let Numerical_integration = Integrate (fun x -> (x*x)*(0.1995*exp((-((x-1.0)*(x-1.0)))/8.0))) -10.0 10.0 number_of_draws;;

(* ---------------------------------- *)
(* mean value theorem for integration *)
(* ---------------------------------- *)
// A normally distributed random generator
let normd = new Normal(0.0, 1.0)

let theta_post = seq { for i in 1 .. number_of_draws do yield theta_mean + (sqrt(theta_var)* normd.Sample()) }

// Lets estimate the expected value of the function theta^2 using MC integration
let theta_sq = theta_post |> Seq.map (fun x -> x * x)

// Seq.length theta_sq |> printfn "Length: %d"
let MC_integration = theta_sq.Mean()

printf "Numerical integration (Trapezoidal's Rule)"
Numerical_integration;;
printf "     "
printf "Monte Carlo integration"
MC_integration;;
// printf "Monte Carlo integration (alternative method)"
// theta_sq1;;
```

# Randomize number generator in F# (Bayesian computing in functional language. )

Over the last few days, I study econometric with a book “Gary Koop Bayesian Econometrics”. This book is NOT over theoretical/mathematical orientated. It is important for me, there are many book which are related with “Bayesian Econometrics”, but all of them are over mathematical. This book is an exceptional, Professor Gary provided all source code on his web-site. He have provided some question sheets which take reader through some basic programming skills relating to course materials. Anyone else would know more about Bayesian Econometrics, you should take a look on his web-site. You would find/learn more than you expected.

I wrote a F# versions for his source code. The purpose of this article was to show the algorithms, and to borrow an implementation in F# and the examples may be interesting for those who like bayesian computing in functional language.

Professor Gary Koop conducted a course for various central-bank all over the world.
====================================================================================
Bayesian Methods for Empirical Macroeconomics
A Course Prepared the Global Initiative Program, Bank of Korea Academy
http://personal.strath.ac.uk/gary.koop/BoK_course.html

```(* ---------------------- *)
(* Question 1 Ex1  *)
(* -------------------------- *)
#r "System.Windows.Forms.DataVisualization.dll"
#r @"../packages/MathNet.Numerics.3.2.1/lib/Net40/MathNet.Numerics.dll"

open System
open System.Net
open System.Windows.Forms
open System.Windows.Forms.DataVisualization.Charting
open Microsoft.FSharp.Control.WebExtensions
open MathNet.Numerics.Distributions
open MathNet.Numerics.Statistics

let num_draws = 10000

(* ---------------------- *)
(* Do results for uniform random number generator *)
(* -------------------------- *)
let uniformRandom =
printf "RESULTS - MEAN, STD "
printf "-----------------------------------------------------------------"
printf " "
printf "UNIFORM"
let normd = new Normal(0.0, 1.0)
let tempp = normd.Samples() |> Seq.take num_draws |> Seq.toList
let unifans = [tempp.Mean(), tempp.StandardDeviation()]
let trueans = [0.5, (1.0/sqrt(12.0))]
printf "%f %f\n",unifans, trueans

(* ---------------------- *)
(* Do results for Normal random number generator  *)
(* ---------------------- *)
let standardNormalRandom =
printf "STANDARD NORMAL"
let normd = new Normal(0.0, 1.0)
let tempp = normd.Samples() |> Seq.take num_draws |> Seq.toList
let normans = [tempp.Mean(), tempp.StandardDeviation()]
let trueans = [0.0, 1.0]
printf "%f %f\n",normans, trueans

(* ---------------------- *)
(* Do results for Student-t random number generator  *)
(* ---------------------- *)
let studentTRandom =
printf "STUDENT-T(3)"
let stdT = new StudentT(3.0,(float num_draws),1.0)
let tempp = stdT.Samples() |> Seq.take num_draws |> Seq.toList
let tans = [tempp.Mean(), tempp.StandardDeviation()]
let trueans = [0.0, sqrt(3.0)]
printf "%f %f\n",tans, trueans

(* ---------------------- *)
(* Do results for Beta random number generator *)
(* ---------------------- *)
let BETARandom =
printf "BETA(3,2)"
let beta = new Beta(3.0,2.0)
let tempp = beta.Samples() |> Seq.take num_draws |> Seq.toList
let betaans = [tempp.Mean(), tempp.StandardDeviation()]
let trueans = [(3/(3+2)),( 1 / 5 )]
printf "%f %f\n",betaans, trueans

(* ---------------------- *)
(* Do results for Exponential random number generator *)
(* ---------------------- *)
let exponentialRandom =
printf "Exponential (5)"
let exp = new Exponential(5.0)
let tempp = exp.Samples() |> Seq.take num_draws |> Seq.toList
let expans = [tempp.Mean(), tempp.StandardDeviation()]
let trueans = [5.0, 5.0]
printf "%f %f\n",expans, trueans

(* ---------------------- *)
(* Do results for Chi-Square random number generator *)
(* ---------------------- *)
let chiSquareRandom =
printf "ChiSquared(3)"
let chi2 = new ChiSquared(3.0)
let tempp = chi2.Samples() |> Seq.take num_draws |> Seq.toList
let chi2ans = [tempp.Mean(), tempp.StandardDeviation()]
let trueans = [3.0, sqrt(6.0)]
printf "%f %f\n",chi2ans, trueans

(* ---------------------- *)
(* Do results for Gamma random number generator *)
(* ---------------------- *)
let gammaRandom =
printf "Gamma(4,2)"
let gamma = new Gamma(4.0,2.0)
let tempp = gamma.Samples() |> Seq.take num_draws |> Seq.toList
let gamans = [tempp.Mean(), tempp.StandardDeviation()]
let trueans = [8.0, 4.0]
printf "%f %f\n",gamans, trueans

uniformRandom;;
standardNormalRandom;;
studentTRandom;;
BETARandom;;
exponentialRandom;;
chiSquareRandom;;
gammaRandom;;

```

# Generate Random number (Brownian motion), Computational Expression F#

To implement a standard Brownian Motion, computation expressions in F# provided a convenient syntax for writing a random number generator; it can be used to manage data, sequence, and expression in a single line of code.

A normally distributed random generator

```
let t = seq { 0.0 .. dt .. T };;

/// Generates random walk starting from value 'v'
let rec randomWalk v =
seq { // Emit the first value of the walk
yield v
// Emit all values starting from new location
yield! randomWalk (v + sqrt(dt) * normd.Sample()) }

// Initial call to generate random walk from value 0
let W = randomWalk(0.0) |> Seq.take (int N) |> Seq.toList

do (Seq.take (int N) W ) |> Seq.iter ( wienerProcess.Points.Add >> ignore)

```