Logo ROOT  
Reference Guide
portfolio.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_quadp
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}
14 /// gSystem->Load("libQuadp");
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 ///
26 /// You might wonder what is so special about this objective which is quadratic in
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 ///
44 /// Enough said about quadratic programming, let's return to our example .
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
110 /// (Google for "slippage optimization")
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 
139 const Int_t nrStocks = 10;
140 static const Char_t *stocks[] =
141  {"GE","SUNW","QCOM","BRCM","TYC","IBM","AMAT","C","PFE","HD"};
142 
143 class TStockDaily {
144 public:
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;
151  Int_t fCloseAdj; // 100*close_price adjusted for splits and dividend
152 
153  TStockDaily() {
154  fDate = fVol = fOpen = fHigh = fLow = fClose = fCloseAdj = 0;
155  }
156  virtual ~TStockDaily() {}
157 
158  ClassDef(TStockDaily,1)
159 };
160 
161 //---------------------------------------------------------------------------
162 Double_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 //---------------------------------------------------------------------------
168 TArrayF &StockReturn(TFile *f,const TString &name,Int_t sDay,Int_t eDay)
169 {
170  TTree *tDaily = (TTree*)f->Get(name);
171  TStockDaily *data = 0;
172  tDaily->SetBranchAddress("daily",&data);
173  TBranch *b_closeAdj = tDaily->GetBranch("fCloseAdj");
174  TBranch *b_date = tDaily->GetBranch("fDate");
175 
176  //read only the "adjusted close" branch for all entries
177  const Int_t nrEntries = (Int_t)tDaily->GetEntries();
178  TArrayF closeAdj(nrEntries);
179  for (Int_t i = 0; i < nrEntries; i++) {
180  b_date->GetEntry(i);
181  b_closeAdj->GetEntry(i);
182  if (data->fDate >= sDay && data->fDate <= eDay)
183  closeAdj[i] = data->fCloseAdj/100.;
184  }
185 
186  TArrayF *r = new TArrayF(nrEntries-1);
187  for (Int_t i = 1; i < nrEntries; i++)
188 // (*r)[i-1] = closeAdj[i]-closeAdj[i-1];
189  (*r)[i-1] = closeAdj[i]/closeAdj[i-1];
190 
191  return *r;
192 }
193 
194 #ifndef __MAKECINT__
195 //---------------------------------------------------------------------------
196 TVectorD 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 
258  TQpDataDens *prob = (TQpDataDens *)qp->MakeData(c,Q,xlo,ixlo,xup,ixup,A,b,C,clo,iclo,cup,icup);
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  //---------------------------------------------------------------------------
285 void 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);
339  gPad->SetGridx();
340  gPad->SetGridy();
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);
356  legend1->AddEntry(f1,"1-exp(-2.0*x)","l");
357  legend1->AddEntry(f2,"1-exp(-10.*x)","l");
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);
381  legend2->AddEntry(h1,"high risk","f");
382  legend2->AddEntry(h2,"low risk","f");
383  legend2->Draw();
384 }
c
#define c(i)
Definition: RSha256.hxx:119
TF1::GetHistogram
virtual TH1 * GetHistogram() const
Return a pointer to the histogram used to visualise the function.
Definition: TF1.cxx:1585
TVectorD.h
TQpProbDens::MakeData
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
f
#define f(i)
Definition: RSha256.hxx:122
TH1::SetMinimum
virtual void SetMinimum(Double_t minimum=-1111)
Definition: TH1.h:396
TQpProbDens::MakeVariables
virtual TQpVar * MakeVariables(const TQpDataBase *data)
Setup the variables.
Definition: TQpProbDens.cxx:182
Form
char * Form(const char *fmt,...)
r
ROOT::R::TRInterface & r
Definition: Object.C:4
TLegend.h
TTree
Definition: TTree.h:79
TMath::Sqrt
Double_t Sqrt(Double_t x)
Definition: TMath.h:680
TTree::SetBranchAddress
virtual Int_t SetBranchAddress(const char *bname, void *add, TBranch **ptr=0)
Change branch address, dealing with clone trees properly.
Definition: TTree.cxx:8205
TMath::Exp
Double_t Exp(Double_t x)
Definition: TMath.h:716
TROOT::GetTutorialDir
static const TString & GetTutorialDir()
Get the tutorials directory in the installation. Static utility function.
Definition: TROOT.cxx:3019
TGeant4Unit::s
static constexpr double s
Definition: TGeant4SystemOfUnits.h:168
TFile::Open
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:3946
Int_t
int Int_t
Definition: RtypesCore.h:45
x
Double_t x[n]
Definition: legend1.C:17
TMatrixTSym
Definition: TMatrixDSymfwd.h:22
ROOT::Math::Cephes::A
static double A[]
Definition: SpecFuncCephes.cxx:170
TGondzioSolver.h
TCanvas.h
TTree.h
TH1::SetXTitle
virtual void SetXTitle(const char *title)
Definition: TH1.h:410
TString
Definition: TString.h:136
TSystem::AccessPathName
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
TArrayF
Definition: TArrayF.h:27
TMatrixT
Definition: TMatrixDfwd.h:22
TH1::SetYTitle
virtual void SetYTitle(const char *title)
Definition: TH1.h:411
b
#define b(i)
Definition: RSha256.hxx:118
TFile.h
h1
TH1F * h1
Definition: legend1.C:5
TQpProbDens::MakeResiduals
virtual TQpResidual * MakeResiduals(const TQpDataBase *data)
Setup the residuals.
Definition: TQpProbDens.cxx:172
TMatrixDSym.h
TF1::SetParameter
virtual void SetParameter(Int_t param, Double_t value)
Definition: TF1.h:630
ROOT::Math::Cephes::C
static double C[]
Definition: SpecFuncCephes.cxx:187
TBranch
Definition: TBranch.h:89
TTree::GetBranch
virtual TBranch * GetBranch(const char *name)
Return pointer to the branch with the given name in this tree or its friends.
Definition: TTree.cxx:5210
TSystem.h
TQpProbDens
Definition: TQpProbDens.h:60
TH1::SetMaximum
virtual void SetMaximum(Double_t maximum=-1111)
Definition: TH1.h:395
TLegend::AddEntry
TLegendEntry * AddEntry(const TObject *obj, const char *label="", Option_t *option="lpf")
Add a new entry to this legend.
Definition: TLegend.cxx:330
TBranch::GetEntry
virtual Int_t GetEntry(Long64_t entry=0, Int_t getall=0)
Read all leaves of entry and return total number of bytes read.
Definition: TBranch.cxx:1582
TH1::Fill
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3274
TQpResidual
Definition: TQpResidual.h:61
TArrayF.h
TH1::SetStats
virtual void SetStats(Bool_t stats=kTRUE)
Set statistics option on/off.
Definition: TH1.cxx:8445
TFile
Definition: TFile.h:54
TQpVar
Definition: TQpVar.h:59
gSystem
R__EXTERN TSystem * gSystem
Definition: TSystem.h:559
sum
static long int sum(long int i)
Definition: Factory.cxx:2272
TVectorT
Definition: TMatrixTBase.h:78
f1
TF1 * f1
Definition: legend1.C:11
Double_t
double Double_t
Definition: RtypesCore.h:59
TGondzioSolver
Definition: TGondzioSolver.h:56
TQpDataDens
Definition: TQpDataDens.h:62
TF1.h
TCanvas
Definition: TCanvas.h:23
TVectorD
TVectorT< Double_t > TVectorD
Definition: TVectorDfwd.h:22
TH1F
1-D histogram with a float per channel (see TH1 documentation)}
Definition: TH1.h:572
TF1::Draw
virtual void Draw(Option_t *option="")
Draw this function with its current attributes.
Definition: TF1.cxx:1338
ClassDef
#define ClassDef(name, id)
Definition: Rtypes.h:325
name
char name[80]
Definition: TGX11.cxx:110
gPad
#define gPad
Definition: TVirtualPad.h:287
TLegend::Draw
virtual void Draw(Option_t *option="")
Draw this legend with its current attributes.
Definition: TLegend.cxx:423
TLegend
Definition: TLegend.h:23
TF1
1-Dim function class
Definition: TF1.h:212
TMatrixD.h
Riostream.h
Char_t
char Char_t
Definition: RtypesCore.h:33
TH1.h
TTree::GetEntries
virtual Long64_t GetEntries() const
Definition: TTree.h:458
TQpVar::fX
TVectorD fX
Definition: TQpVar.h:126
TH1::SetBarWidth
virtual void SetBarWidth(Float_t width=0.5)
Definition: TH1.h:357
TH1::SetBarOffset
virtual void SetBarOffset(Float_t offset=0.25)
Definition: TH1.h:356
TArray::GetSize
Int_t GetSize() const
Definition: TArray.h:47
TMath.h
TQpProbDens.h
TH1::Draw
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2997
c1
return c1
Definition: legend1.C:41
ROOT::Math::Cephes::Q
static double Q[]
Definition: SpecFuncCephes.cxx:294