ROOT   Reference Guide
portfolio.C
Go to the documentation of this file.
1/// \file
3/// \notebook
4/// This macro shows in detail the use of the quadratic programming package quadp .
5/// Running this macro :
6///
7/// ~~~{.cpp}
8/// .x portfolio.C+
9/// ~~~
10///
11/// or
12///
13/// ~~~{.cpp}
15/// .L portFolio.C+; portfolio()
16/// ~~~
17///
18/// Let's first review what we exactly mean by "quadratic programming" :
19///
20/// We want to minimize the following objective function :
21///
22/// \f$c^T x + ( 1/2 ) x^T Q x \f$ wrt. the vector \f$x \f$
23///
24/// \f$c \f$ is a vector and \f$Q \f$ a symmetric positive definite matrix
25///
27/// the unknowns, that can not be done by Minuit/Fumili . Well, we have in addition
28/// the following boundary conditions on \f$x \f$:
29///
30/// \f[
31/// A x = b \\
32/// clo \le C x \le cup \\
33/// xlo \le x \le xup
34/// \f] where A and C are arbitrary matrices and the rest are vectors
35///
36/// Not all these constraints have to be defined . Our example will only use \f$xlo \f$,
37/// \f$A \f$ and \f$b \f$
38/// Still, this could be handled by a general non-linear minimizer like Minuit by introducing
39/// so-called "slack" variables . However, quadp is tailored to objective functions not more
40/// complex than being quadratic . This allows usage of solving techniques which are even
41/// stable for problems involving for instance 500 variables, 100 inequality conditions
42/// and 50 equality conditions .
43///
45/// Suppose, after a long day of doing physics, you have a look at your investments and
46/// realize that an early retirement is not possible, given the returns of your stocks .
47/// So what now ? ROOT to the rescue ...
48///
49/// In 1990 Harry Markowitz was awarded the Nobel prize for economics: " his work provided new tools
50/// for weighing the risks and rewards of different investments and for valuing corporate stocks and bonds" .
51/// In plain English, he developed the tools to balance greed and fear, we want the maximum return
52/// with the minimum amount of risk. Our stock portfolio should be at the
53/// ["Efficient Frontier"](see http://www.riskglossary.com/articles/efficient_frontier.htm).
54/// To quantify better the risk we are willing to take, we define a utility function \f$U(x) \f$. It describes
55/// as a function of our total assets \f$x \f$, our "satisfaction" . A common choice is \f$1-exp(-k*x) \f$ (the reason for
56/// the exponent will be clear later) . The parameter \f$k \f$ is the risk-aversion factor . For small values of \f$k \f$
57/// the satisfaction is small for small values of \f$x \f$; by increasing \f$x \f$ the satisfaction can still be increased
58/// significantly . For large values of \f$k \f$, \f$U(x) \f$ increases rapidly to 1, there is no increase in satisfaction
59/// for additional dollars earned .
60///
61/// In summary :
62/// - small \f$k \f$ ==> risk-loving investor
63/// - large \f$k \f$ ==> risk-averse investor
64///
65/// Suppose we have for nrStocks the historical daily returns \f$r = closing_price(n) - closing_price(n-1) \f$.
66/// Define a vector \f$x \f$ of length of \f$nrStocks \f$, which contains the fraction of our money invested in
67/// each stock . We can calculate the average daily return \f$z \f$ of our portfolio and its variance using
68/// the portfolio covariance Covar :
69///
70/// \f$z = r^T x \f$ and \f$var = x^T Covar x \f$
71///
72/// Assuming that the daily returns have a Normal distribution, \f$N(x) \f$, so will \f$z \f$ with mean \f$r^T x \f$
73/// and variance \f$x^T Covar x \f$
74///
75/// The expected value of the utility function is :
76///
77/// \f[
78/// E(u(x)) = Int (1-exp(-k*x) N(x) dx \\
79/// = 1-exp(-k (r^T x - 0.5 k x^T Covar x) ) \\
80/// \f]
81///
82/// Its value is maximised by maximising \f$r^T x -0.5 k x^T Covar x \f$
83/// under the condition \f$sum (x_i) = 1 \f$, meaning we want all our money invested and
84/// \f$x_i \ge 0 \f$, we can not "short" a stock
85///
86/// For 10 stocks we got the historical daily data for Sep-2000 to Jun-2004:
87///
88/// - GE : General Electric Co
89/// - SUNW : Sun Microsystems Inc
90/// - QCOM : Qualcomm Inc
91/// - BRCM : Broadcom Corp
92/// - TYC : Tyco International Ltd
93/// - IBM : International Business Machines Corp
94/// - AMAT : Applied Materials Inc
95/// - C : Citigroup Inc
96/// - PFE : Pfizer Inc
97/// - HD : Home Depot Inc
98///
99/// We calculate the optimal portfolio for 2.0 and 10.0 .
100///
101/// Food for thought :
102/// - We assumed that the stock returns have a Normal distribution . Check this assumption by
103/// histogramming the stock returns !
104/// - We used for the expected return in the objective function, the flat average over a time
105/// period . Investment firms will put significant resources in improving the return prediction .
106/// - If you want to trade significant number of shares, several other considerations have
107/// to be taken into account :
108/// + If you are going to buy, you will drive the price up (so-called "slippage") .
109/// This can be taken into account by adding terms to the objective
111/// + FTC regulations might have to be added to the inequality constraints
112/// - Investment firms do not want to be exposed to the "market" as defined by a broad
113/// index like the S&P and "hedge" this exposure away . A perfect hedge this can be added
114/// as an equality constrain, otherwise add an inequality constrain .
115///
116/// \macro_image
117/// \macro_output
118/// \macro_code
119///
120/// \author Eddy Offermann
121
122#include "Riostream.h"
123#include "TCanvas.h"
124#include "TFile.h"
125#include "TMath.h"
126#include "TTree.h"
127#include "TArrayF.h"
128#include "TH1.h"
129#include "TF1.h"
130#include "TLegend.h"
131#include "TSystem.h"
132
133#include "TMatrixD.h"
134#include "TMatrixDSym.h"
135#include "TVectorD.h"
136#include "TQpProbDens.h"
137#include "TGondzioSolver.h"
138
139const Int_t nrStocks = 10;
140static const Char_t *stocks[] =
141 {"GE","SUNW","QCOM","BRCM","TYC","IBM","AMAT","C","PFE","HD"};
142
143class TStockDaily {
144public:
145 Int_t fDate;
146 Int_t fOpen; // 100*open_price
147 Int_t fHigh; // 100*high_price
148 Int_t fLow; // 100*low_price
149 Int_t fClose; // 100*close_price
150 Int_t fVol;
152
153 TStockDaily() {
154 fDate = fVol = fOpen = fHigh = fLow = fClose = fCloseAdj = 0;
155 }
156 virtual ~TStockDaily() {}
157
158 ClassDef(TStockDaily,1)
159};
160
161//---------------------------------------------------------------------------
162Double_t RiskProfile(Double_t *x, Double_t *par) {
163 Double_t riskFactor = par[0];
164 return 1-TMath::Exp(-riskFactor*x[0]);
165}
166
167//---------------------------------------------------------------------------
168TArrayF &StockReturn(TFile *f,const TString &name,Int_t sDay,Int_t eDay)
169{
170 TTree *tDaily = (TTree*)f->Get(name);
171 TStockDaily *data = 0;
174 TBranch *b_date = tDaily->GetBranch("fDate");
175
177 const Int_t nrEntries = (Int_t)tDaily->GetEntries();
179 for (Int_t i = 0; i < nrEntries; i++) {
180 b_date->GetEntry(i);
182 if (data->fDate >= sDay && data->fDate <= eDay)
184 }
185
186 TArrayF *r = new TArrayF(nrEntries-1);
187 for (Int_t i = 1; i < nrEntries; i++)
190
191 return *r;
192}
193
194#ifndef __MAKECINT__
195//---------------------------------------------------------------------------
196TVectorD OptimalInvest(Double_t riskFactor,TVectorD r,TMatrixDSym Covar)
197{
198// what the quadratic programming package will do:
199//
200// minimize c^T x + ( 1/2 ) x^T Q x
201// subject to A x = b
202// clo <= C x <= cup
203// xlo <= x <= xup
204// what we want :
205//
206// maximize c^T x - k ( 1/2 ) x^T Q x
207// subject to sum_x x_i = 1
208// 0 <= x_i
209
210 // We have nrStocks weights to determine,
211 // 1 equality- and 0 inequality- equations (the simple square boundary
212 // condition (xlo <= x <= xup) does not count)
213
214 const Int_t nrVar = nrStocks;
215 const Int_t nrEqual = 1;
216 const Int_t nrInEqual = 0;
217
218 // flip the sign of the objective function because we want to maximize
219 TVectorD c = -1.*r;
220 TMatrixDSym Q = riskFactor*Covar;
221
222 // equality equation
223 TMatrixD A(nrEqual,nrVar); A = 1;
224 TVectorD b(nrEqual); b = 1;
225
226 // inequality equation
227 //
228 // - although not applicable in the current situation since nrInEqual = 0, one
229 // has to specify not only clo and cup but also an index vector iclo and icup,
230 // whose values are either 0 or 1 . If iclo[j] = 1, the lower boundary condition
231 // is active on x[j], etc. ...
232
233 TMatrixD C (nrInEqual,nrVar);
234 TVectorD clo (nrInEqual);
235 TVectorD cup (nrInEqual);
236 TVectorD iclo(nrInEqual);
237 TVectorD icup(nrInEqual);
238
239 // simple square boundary condition : 0 <= x_i, so only xlo is relevant .
240 // Like for clo and cup above, we have to define an index vector ixlo and ixup .
241 // Since each variable has the lower boundary, we can set the whole vector
242 // ixlo = 1
243
244 TVectorD xlo (nrVar); xlo = 0;
245 TVectorD xup (nrVar); xup = 0;
246 TVectorD ixlo(nrVar); ixlo = 1;
247 TVectorD ixup(nrVar); ixup = 0;
248
249 // setup the quadratic programming problem . Since a small number of variables are
250 // involved and "Q" has everywhere entries, we chose the dense version "TQpProbDens" .
251 // In case of a sparse formulation, simply replace all "Dens" by "Sparse" below and
252 // use TMatrixDSparse instead of TMatrixDSym and TMatrixD
253
254 TQpProbDens *qp = new TQpProbDens(nrVar,nrEqual,nrInEqual);
255
256 // stuff all the matrices/vectors defined above in the proper places
257
259
260 // setup the nrStock variables, vars->fX will contain the final solution
261
262 TQpVar *vars = qp->MakeVariables(prob);
263 TQpResidual *resid = qp->MakeResiduals(prob);
264
265 // Now we have to choose the method of solving, either TGondzioSolver or TMehrotraSolver
266 // The Gondzio method is more sophisticated and therefore numerically more involved
267 // If one want the Mehrotra method, simply replace "Gondzio" by "Mehrotra" .
268
269 TGondzioSolver *s = new TGondzioSolver(qp,prob);
270 const Int_t status = s->Solve(prob,vars,resid);
271
272 const TVectorD weight = vars->fX;
273
274 delete qp; delete prob; delete vars; delete resid; delete s;
275 if (status != 0) {
276 cout << "Could not solve this problem." <<endl;
277 return TVectorD(nrStocks);
278 }
279
280 return weight;
281}
282#endif
283
284 //---------------------------------------------------------------------------
285void portfolio()
286{
287 const Int_t sDay = 20000809;
288 const Int_t eDay = 20040602;
289
290 const char *fname = "stock.root";
291 TFile *f = 0;
292 if (!gSystem->AccessPathName(fname)) {
293 f = TFile::Open(fname);
294 } else if (!gSystem->AccessPathName(Form("%s/quadp/%s", TROOT::GetTutorialDir().Data(), fname))) {
295 f = TFile::Open(Form("%s/quadp/%s", TROOT::GetTutorialDir().Data(), fname));
296 } else {
297 printf("accessing %s file from http://root.cern.ch/files\n",fname);
298 f = TFile::Open(Form("http://root.cern.ch/files/%s",fname));
299 }
300 if (!f) return;
301
302 TArrayF *data = new TArrayF[nrStocks];
303 for (Int_t i = 0; i < nrStocks; i++) {
304 const TString symbol = stocks[i];
305 data[i] = StockReturn(f,symbol,sDay,eDay);
306 }
307
308 const Int_t nrData = data[0].GetSize();
309
310 TVectorD r(nrStocks);
311 for (Int_t i = 0; i < nrStocks; i++)
312 r[i] = data[i].GetSum()/nrData;
313
314 TMatrixDSym Covar(nrStocks);
315 for (Int_t i = 0; i < nrStocks; i++) {
316 for (Int_t j = 0; j <= i; j++) {
317 Double_t sum = 0.;
318 for (Int_t k = 0; k < nrData; k++) {
319 sum += (data[i][k] - r[i]) * (data[j][k] - r[j]);
320 }
321 Covar(i,j) = Covar(j,i) = sum/nrData;
322 }
323 }
324
325 const TVectorD weight1 = OptimalInvest(2.0,r,Covar);
326 const TVectorD weight2 = OptimalInvest(10.,r,Covar);
327
328 cout << "stock daily daily w1 w2" <<endl;
329 cout << "symb return sdv " <<endl;
330 for (Int_t i = 0; i < nrStocks; i++)
331 printf("%s\t: %.3f %.3f %.3f %.3f\n",stocks[i],r[i],TMath::Sqrt(Covar[i][i]),weight1[i],weight2[i]);
332
333 TCanvas *c1 = new TCanvas("c1","Portfolio Optimizations",10,10,800,900);
334 c1->Divide(1,2);
335
336 // utility function / risk profile
337
338 c1->cd(1);
341
342 TF1 *f1 = new TF1("f1",RiskProfile,0,2.5,1);
343 f1->SetParameter(0,2.0);
344 f1->SetLineColor(49);
345 f1->Draw("AC");
346 f1->GetHistogram()->SetXTitle("dollar");
347 f1->GetHistogram()->SetYTitle("utility");
348 f1->GetHistogram()->SetMinimum(0.0);
349 f1->GetHistogram()->SetMaximum(1.0);
350 TF1 *f2 = new TF1("f2",RiskProfile,0,2.5,1);
351 f2->SetParameter(0,10.);
352 f2->SetLineColor(50);
353 f2->Draw("CSAME");
354
355 TLegend *legend1 = new TLegend(0.50,0.65,0.70,0.82);
358 legend1->Draw();
359
360 // vertical bar chart of portfolio distribution
361
362 c1->cd(2);
363 TH1F *h1 = new TH1F("h1","Portfolio Distribution",nrStocks,0,0);
364 TH1F *h2 = new TH1F("h2","Portfolio Distribution",nrStocks,0,0);
365 h1->SetStats(0);
366 h1->SetFillColor(49);
367 h2->SetFillColor(50);
368 h1->SetBarWidth(0.45);
369 h1->SetBarOffset(0.1);
370 h2->SetBarWidth(0.4);
371 h2->SetBarOffset(0.55);
372 for (Int_t i = 0; i < nrStocks; i++) {
373 h1->Fill(stocks[i],weight1[i]);
374 h2->Fill(stocks[i],weight2[i]);
375 }
376
377 h1->Draw("BAR2 HIST");
378 h2->Draw("BAR2SAME HIST");
379
380 TLegend *legend2 = new TLegend(0.50,0.65,0.70,0.82);
383 legend2->Draw();
384}
ROOT::R::TRInterface & r
Definition: Object.C:4
#define b(i)
Definition: RSha256.hxx:100
#define f(i)
Definition: RSha256.hxx:104
#define c(i)
Definition: RSha256.hxx:101
int Int_t
Definition: RtypesCore.h:45
char Char_t
Definition: RtypesCore.h:37
double Double_t
Definition: RtypesCore.h:59
#define ClassDef(name, id)
Definition: Rtypes.h:325
char name[80]
Definition: TGX11.cxx:110
char * Form(const char *fmt,...)
R__EXTERN TSystem * gSystem
Definition: TSystem.h:559
TVectorT< Double_t > TVectorD
Definition: TVectorDfwd.h:23
Array of floats (32 bits per element).
Definition: TArrayF.h:27
Int_t GetSize() const
Definition: TArray.h:47
virtual void SetFillColor(Color_t fcolor)
Set the fill area color.
Definition: TAttFill.h:37
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition: TAttLine.h:40
A TTree is a list of TBranches.
Definition: TBranch.h:89
virtual Int_t GetEntry(Long64_t entry=0, Int_t getall=0)
Definition: TBranch.cxx:1652
The Canvas class.
Definition: TCanvas.h:23
1-Dim function class
Definition: TF1.h:213
virtual TH1 * GetHistogram() const
Return a pointer to the histogram used to visualise the function Note that this histogram is managed ...
Definition: TF1.cxx:1574
virtual void Draw(Option_t *option="")
Draw this function with its current attributes.
Definition: TF1.cxx:1322
virtual void SetParameter(Int_t param, Double_t value)
Definition: TF1.h:634
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format.
Definition: TFile.h:54
static TFile * Open(const char *name, Option_t *option="", const char *ftitle="", Int_t compress=ROOT::RCompressionSetting::EDefaults::kUseCompiledDefault, Int_t netopt=0)
Create / open a file.
Definition: TFile.cxx:3997
Derived class of TQpSolverBase implementing Gondzio-correction version of Mehrotra's original predict...
1-D histogram with a float per channel (see TH1 documentation)}
Definition: TH1.h:575
virtual void SetBarOffset(Float_t offset=0.25)
Set the bar offset as fraction of the bin width for drawing mode "B".
Definition: TH1.h:359
virtual void SetXTitle(const char *title)
Definition: TH1.h:413
virtual void SetMaximum(Double_t maximum=-1111)
Definition: TH1.h:398
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3350
virtual void SetMinimum(Double_t minimum=-1111)
Definition: TH1.h:399
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:3073
virtual void SetYTitle(const char *title)
Definition: TH1.h:414
virtual void SetBarWidth(Float_t width=0.5)
Set the width of bars as fraction of the bin width for drawing mode "B".
Definition: TH1.h:360
virtual void SetStats(Bool_t stats=kTRUE)
Set statistics option on/off.
Definition: TH1.cxx:8830
This class displays a legend box (TPaveText) containing several legend entries.
Definition: TLegend.h:23
TLegendEntry * AddEntry(const TObject *obj, const char *label="", Option_t *option="lpf")
Add a new entry to this legend.
Definition: TLegend.cxx:330
virtual void Draw(Option_t *option="")
Draw this legend with its current attributes.
Definition: TLegend.cxx:423
Data for the dense QP formulation.
dense matrix problem formulation
Definition: TQpProbDens.h:61
virtual TQpDataBase * MakeData(Double_t *c, Double_t *Q, Double_t *xlo, Bool_t *ixlo, Double_t *xup, Bool_t *ixup, Double_t *A, Double_t *bA, Double_t *C, Double_t *clo, Bool_t *iclo, Double_t *cup, Bool_t *icup)
Setup the data.
Definition: TQpProbDens.cxx:80
virtual TQpResidual * MakeResiduals(const TQpDataBase *data)
Setup the residuals.
virtual TQpVar * MakeVariables(const TQpDataBase *data)
Setup the variables.
The Residuals class calculates and stores the quantities that appear on the right-hand side of the li...
Definition: TQpResidual.h:62
Class containing the variables for the general QP formulation.
Definition: TQpVar.h:60
TVectorD fX
Definition: TQpVar.h:91
static const TString & GetTutorialDir()
Get the tutorials directory in the installation. Static utility function.
Definition: TROOT.cxx:3027
Basic string class.
Definition: TString.h:136
virtual Bool_t AccessPathName(const char *path, EAccessMode mode=kFileExists)
Returns FALSE if one can access a file using the specified access mode.
Definition: TSystem.cxx:1294
A TTree represents a columnar dataset.
Definition: TTree.h:79
virtual TBranch * GetBranch(const char *name)
Return pointer to the branch with the given name in this tree or its friends.
Definition: TTree.cxx:5252
Change branch address, dealing with clone trees properly.
Definition: TTree.cxx:8332
virtual Long64_t GetEntries() const
Definition: TTree.h:459
return c1
Definition: legend1.C:41
Double_t x[n]
Definition: legend1.C:17
TH1F * h1
Definition: legend1.C:5
TF1 * f1
Definition: legend1.C:11
static double Q[]
static double A[]
static double C[]
static constexpr double s
Double_t Exp(Double_t x)
Definition: TMath.h:727
Double_t Sqrt(Double_t x)
Definition: TMath.h:691
static uint64_t sum(uint64_t i)
Definition: Factory.cxx:2345