#include "Riostream.h"
#include "TCanvas.h"
#include "TFile.h"
#include "TMath.h"
#include "TTree.h"
#include "TArrayF.h"
#include "TH1.h"
#include "TF1.h"
#include "TLegend.h"
#include "TSystem.h"

#include "TMatrixD.h"
#include "TMatrixDSym.h"
#include "TVectorD.h"
#include "TQpProbDens.h"
#include "TGondzioSolver.h"

//+ This macro shows in detail the use of the quadratic programming package quadp .
// Running this macro :
//     .x portfolio.C+
//
// Let's first review what we exactly mean by "quadratic programming" :
//
// We want to minimize the following objective function :
//
//   c^T x + ( 1/2 ) x^T Q x    wrt. the vector x
//
//   c is a vector and Q a symmetric positive definite matrix
//
// the unknowns, that can not be done by Minuit/Fumili . Well, we have in addition
// the following boundary conditions on x:
//
//          A x =  b
//  clo <=  C x <= cup
//  xlo <=    x <= xup  ,  where A and C are arbitray matrices and the rest are vectors
//
// Not all these constraints have to be defined . Our example will only use xlo, A and b
// Still, this could be handled by a general non-linear minimizer like Minuit by introducing
// so-called "slack" variables . However, quadp is tailored to objective functions not more
// complex than being quadratic . This allows usage of solving techniques which are even
// stable for problems involving for instance 500 variables, 100 inequality conditions
// and 50 equality conditions .
//
// Suppose, after a long day of doing physics, you have a look at your investments and
// realize that an early retirement is not possible, given the returns of your stocks .
// So what now ? ROOT to the rescue ......
//
// In 1990 Harry Markowitz was awarded the Noble prize for economics : " his work provided new tools
// for weighing the risks and rewards of different investments and for valuing corporate stocks and bonds" .
// In plain English, he developed the tools to balance greed and fear, we want the maximum return
// with the minimum amount of risk. Our stock portfolio should be at the "Efficient Frontier",
// see http://www.riskglossary.com/articles/efficient_frontier.htm .
// To quantify better the risk we are willing to take, we define a utility function U(x) . It describes
// as a function of our total assets x, our "satisfaction" . A common choice is 1-exp(-k*x) (the reason for
// the exponent will be clear later)  . The parameter k is the risk-aversion factor . For small values of k
// the satisfaction is small for small values of x; by increasing x the satisfaction can still be increased
// significantly . For large values of k, U(x) increases rapidly to 1, there is no increase in satisfaction
// for additional dollars earned .
//
// In summary : small k ==> risk-loving investor
//              large k ==> risk-averse investor
//
// Suppose we have for nrStocks the historical daily returns r = closing_price(n) - closing_price(n-1) .
// Define a vector x of length of nrStocks, which contains the fraction of our money invested in
// each stock . We can calculate the average daily return z of our portfolio and its variance using
// the portfolio covariance Covar :
//
//  z = r^T x   and var = x^T Covar x
//
// Assuming that the daily returns have a Normal distribution, N(x), so will z with mean r^T x
// and variance x^T Covar x
//
// The expected value of the utility function is : E(u(x)) = Int (1-exp(-k*x) N(x) dx
//                                                         = 1-exp(-k (r^T x - 0.5 k x^T Covar x) )
//
// Its value is maximized by maximizing  r^T x -0.5 k x^T Covar x
// under the condition sum (x_i) = 1, meaning we want all our money invested and
// x_i >= 0 , we can not "short" a stock
//
// For 10 stocks we got the historical daily data for Sep-2000 to Jun-2004:
//
// GE   : General Electric Co
// SUNW : Sun Microsystems Inc
// QCOM : Qualcomm Inc
// TYC  : Tyco International Ltd
// IBM  : International Business Machines Corp
// AMAT : Applied Materials Inc
// C    : Citigroup Inc
// PFE  : Pfizer Inc
// HD   : Home Depot Inc
//
// We calculate the optimal portfolio for 2.0 and 10.0 .
//
// Food for thought :
// - We assumed that the stock returns have a Normal distribution . Check this assumption by
//   histogramming the stock returns !
// - We used for the expected return in the objective function, the flat average over a time
//   period . Investment firms will put significant resources in improving the return predicton .
// - If you want to trade significant number of shares, several other considerations have
//   to be taken into account :
//    +  If you are going to buy, you will drive the price up (so-called "slippage") .
//       This can be taken into account by adding terms to the objective
//    +  FTC regulations might have to be added to the inequality constraints
// - Investment firms do not want to be exposed to the "market" as defined by a broad
//   index like the S&P and "hedge" this exposure away . A perfect hedge this can be added
//    as an equality constrain, otherwise add an inequality constrain .
//
//Author: Eddy Offermann

const Int_t nrStocks = 10;
static const Char_t *stocks[] =
{"GE","SUNW","QCOM","BRCM","TYC","IBM","AMAT","C","PFE","HD"};

class TStockDaily  {
public:
Int_t fDate;
Int_t fOpen;     // 100*open_price
Int_t fHigh;     // 100*high_price
Int_t fLow;      // 100*low_price
Int_t fClose;    // 100*close_price
Int_t fVol;

TStockDaily() {
fDate = fVol = fOpen = fHigh = fLow = fClose = fCloseAdj = 0;
}
virtual ~TStockDaily() {}

ClassDef(TStockDaily,1)
};

//---------------------------------------------------------------------------
Double_t RiskProfile(Double_t *x, Double_t *par) {
Double_t riskFactor = par[0];
return 1-TMath::Exp(-riskFactor*x[0]);
}

//---------------------------------------------------------------------------
TArrayF &StockReturn(TFile *f,const TString &name,Int_t sDay,Int_t eDay)
{
TTree *tDaily = (TTree*)f->Get(name);
TStockDaily *data = 0;
TBranch *b_date     = tDaily->GetBranch("fDate");

const Int_t nrEntries = (Int_t)tDaily->GetEntries();
for (Int_t i = 0; i < nrEntries; i++) {
b_date->GetEntry(i);
if (data->fDate >= sDay && data->fDate <= eDay)
#ifdef __CINT__
#else
#endif
}

TArrayF *r = new TArrayF(nrEntries-1);
for (Int_t i = 1; i < nrEntries; i++)
#ifdef __CINT__
#else
#endif

return *r;
}

#ifndef __MAKECINT__
//---------------------------------------------------------------------------
TVectorD OptimalInvest(Double_t riskFactor,TVectorD r,TMatrixDSym Covar)
{
// what the quadratic programming package will do:
//
//  minimize    c^T x + ( 1/2 ) x^T Q x
//  subject to                A x  = b
//                    clo <=  C x <= cup
//                    xlo <=    x <= xup
// what we want :
//
//  maximize    c^T x - k ( 1/2 ) x^T Q x
//  subject to        sum_x x_i = 1
//                    0 <= x_i

// We have nrStocks weights to determine,
// 1 equality- and 0 inequality- equations (the simple square boundary
// condition (xlo <= x <= xup) does not count)

const Int_t nrVar     = nrStocks;
const Int_t nrEqual   = 1;
const Int_t nrInEqual = 0;

// flip the sign of the objective function because we want to maximize
TVectorD    c = -1.*r;
TMatrixDSym Q = riskFactor*Covar;

// equality equation
TMatrixD A(nrEqual,nrVar); A = 1;
TVectorD b(nrEqual); b = 1;

// inequality equation
//
// - although not applicable in the current situatio since nrInEqual = 0, one
//   has to specify not only clo and cup but also an index vector iclo and icup,
//   whose values are either 0 or 1 . If iclo[j] = 1, the lower boundary condition
//   is active on x[j], etc. ...

TMatrixD C   (nrInEqual,nrVar);
TVectorD clo (nrInEqual);
TVectorD cup (nrInEqual);
TVectorD iclo(nrInEqual);
TVectorD icup(nrInEqual);

// simple square boundary condition : 0 <= x_i, so only xlo is relevant .
// Like for clo and cup above, we have to define an index vector ixlo and ixup .
// Since each variable has the lower boundary, we can set the whole vector
// ixlo = 1

TVectorD xlo (nrVar); xlo  = 0;
TVectorD xup (nrVar); xup  = 0;
TVectorD ixlo(nrVar); ixlo = 1;
TVectorD ixup(nrVar); ixup = 0;

// setup the quadratic programming problem . Since a small number of variables are
// involved and "Q" has everywhere entries, we chose the dense version "TQpProbDens" .
// In case of a sparse formulation, simply replace all "Dens" by "Sparse" below and
// use TMatrixDSparse instead of TMatrixDSym and TMatrixD

TQpProbDens *qp = new TQpProbDens(nrVar,nrEqual,nrInEqual);

// stuff all the matrices/vectors defined above in the proper places

// setup the nrStock variables, vars->fX will contain the final solution

TQpVar      *vars  = qp->MakeVariables(prob);
TQpResidual *resid = qp->MakeResiduals(prob);

// Now we have to choose the method of solving, either TGondzioSolver or TMehrotraSolver
// The Gondzio method is more sophisticated and therefore numerically more involved
// If one want the Mehrotra method, simply replace "Gondzio" by "Mehrotra" .

TGondzioSolver *s = new TGondzioSolver(qp,prob);
const Int_t status = s->Solve(prob,vars,resid);

const TVectorD weight = vars->fX;

delete qp; delete prob; delete vars; delete resid; delete s;
if (status != 0) {
cout << "Could not solve this problem." <<endl;
return TVectorD(nrStocks);
}

return weight;
}
#endif

//---------------------------------------------------------------------------
void portfolio()
{
const Int_t sDay = 20000809;
const Int_t eDay = 20040602;

const char *fname = "stock.root";
TFile *f = 0;
if (!gSystem->AccessPathName(fname)) {
f = TFile::Open(fname);
} else {
printf("accessing %s file from http://root.cern.ch/files\n",fname);
f = TFile::Open(Form("http://root.cern.ch/files/%s",fname));
}
if (!f) return;

TArrayF *data = new TArrayF[nrStocks];
for (Int_t i = 0; i < nrStocks; i++) {
const TString symbol = stocks[i];
data[i] = StockReturn(f,symbol,sDay,eDay);
}

const Int_t nrData = data[0].GetSize();

TVectorD r(nrStocks);
for (Int_t i = 0; i < nrStocks; i++)
r[i] = data[i].GetSum()/nrData;

TMatrixDSym Covar(nrStocks);
for (Int_t i = 0; i < nrStocks; i++) {
for (Int_t j = 0; j <= i; j++) {
Double_t sum = 0.;
for (Int_t k = 0; k < nrData; k++)
sum += (data[i][k]-r[i])*(data[j][k]-r[j]);
Covar(i,j) = Covar(j,i) = sum/nrData;
}
}

const TVectorD weight1 = OptimalInvest(2.0,r,Covar);
const TVectorD weight2 = OptimalInvest(10.,r,Covar);

cout << "stock     daily  daily   w1     w2" <<endl;
cout << "symb      return sdv              " <<endl;
for (Int_t i = 0; i < nrStocks; i++)
printf("%s\t: %.3f  %.3f  %.3f  %.3f\n",stocks[i],r[i],TMath::Sqrt(Covar[i][i]),weight1[i],weight2[i]);

TCanvas *c1 = new TCanvas("c1","Portfolio Optimizations",10,10,800,900);
c1->Divide(1,2);

// utility function / risk profile

c1->cd(1);

TF1 *f1 = new TF1("f1",RiskProfile,0,2.5,1);
f1->SetParameter(0,2.0);
f1->SetLineColor(49);
f1->Draw("AC");
f1->GetHistogram()->SetXTitle("dollar");
f1->GetHistogram()->SetYTitle("utility");
f1->GetHistogram()->SetMinimum(0.0);
f1->GetHistogram()->SetMaximum(1.0);
TF1 *f2 = new TF1("f2",RiskProfile,0,2.5,1);
f2->SetParameter(0,10.);
f2->SetLineColor(50);
f2->Draw("CSAME");

TLegend *legend1 = new TLegend(0.50,0.65,0.70,0.82);
legend1->Draw();

// vertical bar chart of portfolio distribution

c1->cd(2);
TH1F *h1 = new TH1F("h1","Portfolio Distribution",nrStocks,0,0);
TH1F *h2 = new TH1F("h2","Portfolio Distribution",nrStocks,0,0);
h1->SetStats(0);
h1->SetFillColor(49);
h2->SetFillColor(50);
h1->SetBarWidth(0.45);
h1->SetBarOffset(0.1);
h2->SetBarWidth(0.4);
h2->SetBarOffset(0.55);
for (Int_t i = 0; i < nrStocks; i++) {
h1->Fill(stocks[i],weight1[i]);
h2->Fill(stocks[i],weight2[i]);
}

h1->Draw("BAR2");
h2->Draw("BAR2SAME");

TLegend *legend2 = new TLegend(0.50,0.65,0.70,0.82);
legend2->Draw();
}
portfolio.C:1
portfolio.C:2
portfolio.C:3
portfolio.C:4
portfolio.C:5
portfolio.C:6
portfolio.C:7
portfolio.C:8
portfolio.C:9
portfolio.C:10
portfolio.C:11
portfolio.C:12
portfolio.C:13
portfolio.C:14
portfolio.C:15
portfolio.C:16
portfolio.C:17
portfolio.C:18
portfolio.C:19
portfolio.C:20
portfolio.C:21
portfolio.C:22
portfolio.C:23
portfolio.C:24
portfolio.C:25
portfolio.C:26
portfolio.C:27
portfolio.C:28
portfolio.C:29
portfolio.C:30
portfolio.C:31
portfolio.C:32
portfolio.C:33
portfolio.C:34
portfolio.C:35
portfolio.C:36
portfolio.C:37
portfolio.C:38
portfolio.C:39
portfolio.C:40
portfolio.C:41
portfolio.C:42
portfolio.C:43
portfolio.C:44
portfolio.C:45
portfolio.C:46
portfolio.C:47
portfolio.C:48
portfolio.C:49
portfolio.C:50
portfolio.C:51
portfolio.C:52
portfolio.C:53
portfolio.C:54
portfolio.C:55
portfolio.C:56
portfolio.C:57
portfolio.C:58
portfolio.C:59
portfolio.C:60
portfolio.C:61
portfolio.C:62
portfolio.C:63
portfolio.C:64
portfolio.C:65
portfolio.C:66
portfolio.C:67
portfolio.C:68
portfolio.C:69
portfolio.C:70
portfolio.C:71
portfolio.C:72
portfolio.C:73
portfolio.C:74
portfolio.C:75
portfolio.C:76
portfolio.C:77
portfolio.C:78
portfolio.C:79
portfolio.C:80
portfolio.C:81
portfolio.C:82
portfolio.C:83
portfolio.C:84
portfolio.C:85
portfolio.C:86
portfolio.C:87
portfolio.C:88
portfolio.C:89
portfolio.C:90
portfolio.C:91
portfolio.C:92
portfolio.C:93
portfolio.C:94
portfolio.C:95
portfolio.C:96
portfolio.C:97
portfolio.C:98
portfolio.C:99
portfolio.C:100
portfolio.C:101
portfolio.C:102
portfolio.C:103
portfolio.C:104
portfolio.C:105
portfolio.C:106
portfolio.C:107
portfolio.C:108
portfolio.C:109
portfolio.C:110
portfolio.C:111
portfolio.C:112
portfolio.C:113
portfolio.C:114
portfolio.C:115
portfolio.C:116
portfolio.C:117
portfolio.C:118
portfolio.C:119
portfolio.C:120
portfolio.C:121
portfolio.C:122
portfolio.C:123
portfolio.C:124
portfolio.C:125
portfolio.C:126
portfolio.C:127
portfolio.C:128
portfolio.C:129
portfolio.C:130
portfolio.C:131
portfolio.C:132
portfolio.C:133
portfolio.C:134
portfolio.C:135
portfolio.C:136
portfolio.C:137
portfolio.C:138
portfolio.C:139
portfolio.C:140
portfolio.C:141
portfolio.C:142
portfolio.C:143
portfolio.C:144
portfolio.C:145
portfolio.C:146
portfolio.C:147
portfolio.C:148
portfolio.C:149
portfolio.C:150
portfolio.C:151
portfolio.C:152
portfolio.C:153
portfolio.C:154
portfolio.C:155
portfolio.C:156
portfolio.C:157
portfolio.C:158
portfolio.C:159
portfolio.C:160
portfolio.C:161
portfolio.C:162
portfolio.C:163
portfolio.C:164
portfolio.C:165
portfolio.C:166
portfolio.C:167
portfolio.C:168
portfolio.C:169
portfolio.C:170
portfolio.C:171
portfolio.C:172
portfolio.C:173
portfolio.C:174
portfolio.C:175
portfolio.C:176
portfolio.C:177
portfolio.C:178
portfolio.C:179
portfolio.C:180
portfolio.C:181
portfolio.C:182
portfolio.C:183
portfolio.C:184
portfolio.C:185
portfolio.C:186
portfolio.C:187
portfolio.C:188
portfolio.C:189
portfolio.C:190
portfolio.C:191
portfolio.C:192
portfolio.C:193
portfolio.C:194
portfolio.C:195
portfolio.C:196
portfolio.C:197
portfolio.C:198
portfolio.C:199
portfolio.C:200
portfolio.C:201
portfolio.C:202
portfolio.C:203
portfolio.C:204
portfolio.C:205
portfolio.C:206
portfolio.C:207
portfolio.C:208
portfolio.C:209
portfolio.C:210
portfolio.C:211
portfolio.C:212
portfolio.C:213
portfolio.C:214
portfolio.C:215
portfolio.C:216
portfolio.C:217
portfolio.C:218
portfolio.C:219
portfolio.C:220
portfolio.C:221
portfolio.C:222
portfolio.C:223
portfolio.C:224
portfolio.C:225
portfolio.C:226
portfolio.C:227
portfolio.C:228
portfolio.C:229
portfolio.C:230
portfolio.C:231
portfolio.C:232
portfolio.C:233
portfolio.C:234
portfolio.C:235
portfolio.C:236
portfolio.C:237
portfolio.C:238
portfolio.C:239
portfolio.C:240
portfolio.C:241
portfolio.C:242
portfolio.C:243
portfolio.C:244
portfolio.C:245
portfolio.C:246
portfolio.C:247
portfolio.C:248
portfolio.C:249
portfolio.C:250
portfolio.C:251
portfolio.C:252
portfolio.C:253
portfolio.C:254
portfolio.C:255
portfolio.C:256
portfolio.C:257
portfolio.C:258
portfolio.C:259
portfolio.C:260
portfolio.C:261
portfolio.C:262
portfolio.C:263
portfolio.C:264
portfolio.C:265
portfolio.C:266
portfolio.C:267
portfolio.C:268
portfolio.C:269
portfolio.C:270
portfolio.C:271
portfolio.C:272
portfolio.C:273
portfolio.C:274
portfolio.C:275
portfolio.C:276
portfolio.C:277
portfolio.C:278
portfolio.C:279
portfolio.C:280
portfolio.C:281
portfolio.C:282
portfolio.C:283
portfolio.C:284
portfolio.C:285
portfolio.C:286
portfolio.C:287
portfolio.C:288
portfolio.C:289
portfolio.C:290
portfolio.C:291
portfolio.C:292
portfolio.C:293
portfolio.C:294
portfolio.C:295
portfolio.C:296
portfolio.C:297
portfolio.C:298
portfolio.C:299
portfolio.C:300
portfolio.C:301
portfolio.C:302
portfolio.C:303
portfolio.C:304
portfolio.C:305
portfolio.C:306
portfolio.C:307
portfolio.C:308
portfolio.C:309
portfolio.C:310
portfolio.C:311
portfolio.C:312
portfolio.C:313
portfolio.C:314
portfolio.C:315
portfolio.C:316
portfolio.C:317
portfolio.C:318
portfolio.C:319
portfolio.C:320
portfolio.C:321
portfolio.C:322
portfolio.C:323
portfolio.C:324
portfolio.C:325
portfolio.C:326
portfolio.C:327
portfolio.C:328
portfolio.C:329
portfolio.C:330
portfolio.C:331
portfolio.C:332
portfolio.C:333
portfolio.C:334
portfolio.C:335
portfolio.C:336
portfolio.C:337
portfolio.C:338
portfolio.C:339
portfolio.C:340
portfolio.C:341
portfolio.C:342
portfolio.C:343
portfolio.C:344
portfolio.C:345
portfolio.C:346
portfolio.C:347
portfolio.C:348
portfolio.C:349
portfolio.C:350
portfolio.C:351
portfolio.C:352
portfolio.C:353
portfolio.C:354
portfolio.C:355
portfolio.C:356
portfolio.C:357
portfolio.C:358
portfolio.C:359
portfolio.C:360
portfolio.C:361
portfolio.C:362
portfolio.C:363
portfolio.C:364
portfolio.C:365
portfolio.C:366