Merge pull request #58 from IainNZ/jump010

WIP: JuMPeR v0.2

.travis.yml

@@ -1,23 +1,24 @@
-language: cpp
-compiler:
-    - gcc
+language: julia
+os:
+    - linux
+julia:
+    - release
+    - nightly
 notifications:
     email: false
-env:
-    - JULIAVERSION="juliareleases"
-    - JULIAVERSION="julianightlies"
-before_install:
-    - sudo add-apt-repository ppa:staticfloat/julia-deps -y
-    - sudo add-apt-repository ppa:staticfloat/${JULIAVERSION} -y
-    - sudo apt-get update -qq -y
-    - sudo apt-get install libpcre3-dev libgmp-dev julia -y
-    - git config --global user.name "Travis User"
-    - git config --global user.email "travis@example.net"
-    - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi
+sudo: false
+addons:
+    apt_packages:
+        - gfortran
+        - liblapack-dev
+        - libgmp-dev
+        - libglpk-dev
 script:
-    - julia -e 'versioninfo(); Pkg.init(); Pkg.clone(pwd())'
+    - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi
+    - julia -e 'Pkg.add("JuMP"); Pkg.checkout("JuMP")'
+    - julia -e 'Pkg.clone(pwd())'
     - julia -e 'Pkg.test("JuMPeR", coverage=true)'
-    - if [ $JULIAVERSION = "juliareleases" ]; then julia -e 'Pkg.clone("https://github.com/vgupta1/DDUS.jl.git")'; fi
-    - if [ $JULIAVERSION = "juliareleases" ]; then julia -e 'Pkg.test("DDUS")'; fi
 after_success:
-    - julia -e 'cd(Pkg.dir("JuMPeR")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(Coveralls.process_folder())'
+    - julia -e 'Pkg.add("Coverage")'
+    - julia -e 'cd(Pkg.dir("JuMPeR")); using Coverage; Coveralls.submit(process_folder())'
+    - julia -e 'cd(Pkg.dir("JuMPeR")); using Coverage; Codecov.submit(process_folder())'

README.md

@@ -1,6 +1,6 @@
 ![JuMPeR Logo](http://iainnz.github.io/JuMPeR.jl/logo.svg)
 
-[![Build Status](https://travis-ci.org/IainNZ/JuMPeR.jl.png?branch=master)](https://travis-ci.org/IainNZ/JuMPeR.jl)
+[![Build Status](https://travis-ci.org/IainNZ/JuMPeR.jl.svg?branch=master)](https://travis-ci.org/IainNZ/JuMPeR.jl)
 [![Coverage Status](https://img.shields.io/coveralls/IainNZ/JuMPeR.jl.svg)](https://coveralls.io/r/IainNZ/JuMPeR.jl?branch=master)
 [![JuMPeR](http://pkg.julialang.org/badges/JuMPeR_release.svg)](http://pkg.julialang.org/?pkg=JuMPeR&ver=release)
 

REQUIRE

@@ -1,3 +1,3 @@
 julia 0.3
-JuMP 0.9
+JuMP 0.10 0.11
 Compat

example/portfolio.jl

@@ -1,121 +1,149 @@
-#############################################################################
-# JuMPeR
-# Julia for Mathematical Programming - extension for Robust Optimization
-# See http://github.com/IainNZ/JuMPeR.jl
-#############################################################################
-# portfolio.jl
-# Use the CovarOracle defined in oracle_covar.jl to solve an asset
-# allocation problem. Solves for multiple values of Γ and shows the
-# distribution of returns of the optimal portfolio for each value.
-#############################################################################
-
-using JuMPeR
+#-----------------------------------------------------------------------
+# JuMPeR  --  JuMP Extension for Robust Optimization
+# http://github.com/IainNZ/JuMPeR.jl
+#-----------------------------------------------------------------------
+# Copyright (c) 2015: Iain Dunning
+# This Source Code Form is subject to the terms of the Mozilla Public
+# License, v. 2.0. If a copy of the MPL was not distributed with this
+# file, You can obtain one at http://mozilla.org/MPL/2.0/.
+#-----------------------------------------------------------------------
+# example/portfolio.jl
+# Solve a robust portfolio optimization problem. Can choose between:
+# - Polyhedral or ellipsoidal uncertainty sets
+# - Reformulation or cutting planes
+#-----------------------------------------------------------------------
+
+#-----------------------------------------------------------------------
+# Load packages
+# Modelling
+using JuMP, JuMPeR
+# For generating data
 using Distributions
-using Gurobi
-include("oracle_covar.jl")
+# Only needed for Julia 0.3, used for factorizing covariance matrix
+using Compat
 
+#-----------------------------------------------------------------------
+# Global configuration
+# Number of stocks
 const NUM_ASSET = 10
-
+# Number of samples of past returns
+const NUM_SAMPLES = 1000
+# Uncertainty set type (:Polyhedral or :Ellipsoidal)
+#const UNCERTAINY_SET = :Polyhedral
+const UNCERTAINY_SET = :Ellipsoidal
 # Seed the RNG for reproducibilities sake
-srand(10)
+srand(1234)
 
-#############################################################################
+#-----------------------------------------------------------------------
 # generate_data
-# Sample returns of the assets. The parameters of the assets are fixed, we
-# simply make num_samples draws from the joint distribution.
-# Returns matrix with the samples in rows, each column is an assets
-function generate_data(num_samples::Int)
-    data = zeros(num_samples, NUM_ASSET)
+# Sample returns of the assets from a multivariate normal distribution.
+# Returns a matrix with the samples in rows & each column is an asset.
+function generate_data()
+    data = zeros(NUM_SAMPLES, NUM_ASSET)
 
-    # Linking factors
-    beta = [(i-1.0)/NUM_ASSET for i = 1:NUM_ASSET] 
+    # Linking factors to induce correlations
+    β = [(i-1.0)/NUM_ASSET for i in 1:NUM_ASSET] 
 
-    for sample_ind = 1:num_samples
-        # Common market factor, mean 3%, sd 5%, truncate at +- 3 sd
+    for i in 1:NUM_SAMPLES
+        # Common market factor, μ=3%, σ=5%, truncated at ±3σ
         z = rand(Normal(0.03, 0.05))
         z = max(z, 0.03 - 3*0.05)
         z = min(z, 0.03 + 3*0.05)
 
-        for asset_ind = 1:NUM_ASSET
-            # Idiosyncratic contribution, mean 0%, sd 5%, truncated at +- 3 sd
-            asset = rand(Normal(0.00, 0.05))
-            asset = max(asset, 0.00 - 3*0.05)
-            asset = min(asset, 0.00 + 3*0.05)
-            data[sample_ind, asset_ind] = beta[asset_ind] * z + asset
+        for j in 1:NUM_ASSET
+            # Idiosyncratic contribution, μ=0%, σ=5%, truncated at ±3σ
+            r = rand(Normal(0.00, 0.05))
+            r = max(r, 0.00 - 3*0.05)
+            r = min(r, 0.00 + 3*0.05)
+            data[i,j] = β[j] * z + r
         end
     end
 
     return data
 end
 
-
-#############################################################################
+#-----------------------------------------------------------------------
 # solve_portfolio
-# Solve the robust portfolio problem given a matrix of past returns and
-# a degree of conservatism, Γ. The particular problem solved is
+# Solve the robust portfolio problem given a matrix of past returns.
+# The particular problem solved is
 #   max   z
-#   s.t.  ∑ xi    = 1
-#         ∑ ri*xi ≥ z  ∀ r ∈ R
-#         xi ≥ 0
-#   where R is CovarOracle's set, an ellipse centered on the mean return
-#         and tilted using the covariance of the returns
+#   s.t.  ∑ rᵢ xᵢ ≥ z  ∀ r ∈ R      # Worst-case return
+#         ∑ xᵢ    = 1               # Is a portfolio
+#         xᵢ ≥ 0
+# where R is one of the uncertainty sets:
+# R = { (r,z) | r = μ + Σ^½ z,
+#        Polyhedral:  ‖z‖₁ ≤ Γ, ‖z‖∞ ≤ 1
+#       Ellipsoidal:  ‖z‖₂ ≤ Γ
+#     }
 function solve_portfolio(past_returns, Γ)
     # Setup the robust optimization model
-    m = RobustModel(solver=GurobiSolver(OutputFlag=0))
-
-    # Create the CovarOracle
-    setDefaultOracle!(m, CovarOracle(past_returns, Γ))
-
-    # Variables
-    @defVar(m, obj)  # Put objective as constraint (epigraph form)
-    @defVar(m, 0 <= x[1:NUM_ASSET] <= 1)  # Fractional allocation
+    m = RobustModel()
 
-    # Uncertainties
-    @defUnc(m, r[1:NUM_ASSET])  # The returns
+    # Each asset is a share of the money to invest...
+    @defVar(m, 0 <= x[1:NUM_ASSET] <= 1)
+    # ... and we must allocate all the money.
+    @addConstraint(m, sum(x) == 1)
 
+    # JuMPeR doesn't support uncertain objectives directly.
+    # Uncertain objectives, should put the objective function as
+    # a constraint (epigraph form)
+    @defVar(m, obj)
     @setObjective(m, Max, obj)
 
-    # Portfolio constraint
-    @addConstraint(m, sum(x) == 1)
+    # The uncertain parameters are the returns for each asset
+    @defUnc(m, r[1:NUM_ASSET])
+
+    # The uncertainty set requires the "square root" of the covariance
+    μ = vec(mean(past_returns, 1))  # Want a column vector 
+    Σ = cov(past_returns)
+    L = full(@compat chol(Σ, Val{:L}))  # Σ^½
+    # Define auxiliary uncertain parameters to model underlying factors
+    @defUnc(m, z[1:NUM_ASSET])
+    # Link r with z
+    @addConstraint(m, r .== L*z + μ)
+    if UNCERTAINY_SET == :Polyhedral
+        @addConstraint(m, norm(z,   1) <= Γ)  # ‖z‖₁ ≤ Γ
+        @addConstraint(m, norm(z, Inf) <= 1)  # ‖z‖∞ ≤ 1
+    else
+        @addConstraint(m, norm(z,   2) <= Γ)  # ‖z‖₂ ≤ Γ
+    end
 
-    # The objective constraint - uncertain
-    @addConstraint(m, obj - dot(r, x) <= 0)
+    # The objective function: the worst-case return
+    @addConstraint(m, obj <= dot(r, x))
 
-    # Solve it, report statistics on number of cutting planes etc.
-    #println("Solving model for Γ=$Γ")
-    #solveRobust(m, report=true)
-    #println("Objective value: ", getValue(obj))
-    solveRobust(m)
+    # Solve it
+    solve(m)
 
+    # Return the allocation and the worst-case return
     return getValue(x), getValue(obj)
 end
 
-
-#############################################################################
+#-----------------------------------------------------------------------
 # simulate
-# Generate some past returns, then obtain portfolios for different values
-# of Γ to see the effect of robustness on the distribution of returns
-function simulate(num_past, num_future)
+# Generate some past returns, obtain portfolios for different
+# values of Γ, and analyze distribution of returns
+function simulate()
     # Generate the simulated data
-    past_returns   = generate_data(num_past)
-    future_returns = generate_data(num_future)
+    past_returns   = generate_data()
+    future_returns = generate_data()
 
     # Print table header
-    println("   Γ|   Min|   10%|   20%|  Mean|   Max")
-    for Γ in [0.0, 1.0, 2.0, 3.0, 4.0, 5.0]
+    println("    Γ |    Min |    10% |    20% |   Mean |    Max | Support")
+    for Γ in 0:NUM_ASSET
         # Solve for portfolio
         x, obj = solve_portfolio(past_returns, Γ)
         # Evaluate portfolio returns in future
-        future_z = future_returns*x[:]
+        future_z = future_returns*x
         sort!(future_z)
         min_z    = future_z[1]*100
-        ten_z    = future_z[int(num_future*0.1)]*100
-        twenty_z = future_z[int(num_future*0.2)]*100
+        ten_z    = future_z[div(NUM_SAMPLES,10)]*100
+        twenty_z = future_z[div(NUM_SAMPLES, 5)]*100
         mean_z   = mean(future_z)*100
         max_z    = future_z[end]*100
-        @printf(" %3.1f| %5.1f| %5.1f| %5.1f| %5.1f| %5.1f\n",
-                    Γ, min_z, ten_z, twenty_z, mean_z, max_z)
+        support  = sum(x .> 1e-6)  # Number of assets used
+        @printf(" %4.1f | %6.2f | %6.2f | %6.2f | %6.2f | %6.2f | %7d\n",
+                    Γ, min_z, ten_z, twenty_z, mean_z, max_z, support)
     end
 end
 
-simulate(1000, 1000)
+simulate()