// Black-Litterman metrics example code for SciLab (bl_metrics.sci)
// Copyright (c) Jay Walters, blacklitterman.org, 2009.
//
// Redistribution and use in source and binary forms, 
// with or without modification, are permitted provided 
// that the following conditions are met:
//
// Redistributions of source code must retain the above 
// copyright notice, this list of conditions and the following 
// disclaimer.
// 
// Redistributions in binary form must reproduce the above 
// copyright notice, this list of conditions and the following 
// disclaimer in the documentation and/or other materials 
// provided with the distribution.
//  
// Neither the name of blacklitterman.org nor the names of its
// contributors may be used to endorse or promote products 
// derived from this software without specific prior written
// permission.
//  
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND 
// CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, 
// INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF 
// MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
// DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR 
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 
// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, 
// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR 
// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, 
// WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 
// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH 
// DAMAGE.
//
// This program uses the examples from the paper "The Intuition 
// Behind Black-Litterman Model  Portfolios", by He and Litterman,
// 1999.  You can find a copy of this  paper at the following url.
//     http://papers.ssrn.com/sol3/papers.cfm?abstract_id=334304
//
// For more details on the Black-Litterman model you can also view
// "The BlackLitterman Model: A Detailed Exploration", by this author
// at the following url.
//     http://www.blacklitterman.org/Black-Litterman.pdf
//

// blacklitterman
//   This function performs the Black-Litterman blending of the prior
//   and the views into a new posterior estimate of the returns.
// Inputs
//   delta  - Risk tolerance from the equilibrium portfolio
//   weq    - Weights of the assets in the equilibrium portfolio
//   sigma  - Prior covariance matrix
//   tau    - Coefficiet of uncertainty in the prior estimate of the mean (pi)
//   P      - Pick matrix for the view(s)
//   Q      - Vector of view returns
//   Omega  - Matrix of variance of the views (diagonal)
// Outputs
//   Er     - Posterior estimate of the mean returns
//   pi     - Equilibrium returns used as prior
//   w      - Unconstrained weights computed given the Posterior estimates
//            of the mean and covariance of returns.
//   lambda - A measure of the impact of each view on the posterior estimates.
//
function [er, pi, w, postV, lambda] = blacklitterman(delta, weq, sigma, tau, P, Q, Omega)
  // Reverse optimize and back out the equilibrium returns
  // This is formula (12) page 6.
  pi = weq * sigma * delta;
  // We use tau * sigma many places so just compute it once
  ts = tau * sigma;
  // Compute posterior estimate of the mean
  // This is a simplified version of formula (8) on page 4.
  er = pi' + ts * P' * inv(P * ts * P' + Omega) * (Q - P * pi');
  // Compute posterior estimate of the uncertainty in the mean
  // This is a simplified and combined version of formulas (9) and (15)
  postV = ts - ts * P' * inv(P * ts * P' + Omega) * P * ts;
  posteriorSigma = sigma + postV;
  // Compute posterior weights based on uncertainty in mean
  // May neeed to apply a non-negative constraint to the weights here
  // in order to get it to look pretty.
  w = (er' * inv(delta * posteriorSigma))';
  // We actually will deform the solution with a budget constraint and no
  // short selling to try and present a more realistic example
  C=[ones(pi)]
  ci=zeros(pi)'+0.000001
  // minimizing so p is negative and Q is positive
  //[w,lagr,f]=quapro(delta*sigma,-er,C,[1],ci,[],1);
  // Compute lambda value
  // We solve for lambda from formula (17) page 7, rather than formula (18)
  // just because it is less to type, and we've already computed w*.
  lambda = pinv(P)' * (w'*(1+tau) - weq)';
endfunction

// Function to compute the pdf of the chia square distribution
function [pdf] = pdfchi(x,v)
  pdf = (x ^ (v/2-1)*exp(-x/2))/(2^(v/2)*gamma(v/2));
endfunction

// computeCons
// Computes the Consistency measure from Fusai and Meucci (2003)
//   pi      - Prior estimates of mean
//   Er      - Posterior estimate of mean
//   V       - Covariance of prior estimate of mean
//   P       - Pick matrix for the view(s)
// Outputs
//   md      - Mahalanobis distance of posterior returns from prior returns
//   mdP     - Confidence level that posterior returns are consistent with prior
//   partial - Partial sensitivities of the confidence to the views
//  
function [md,mdP,partial]=computeCons(pi,Er,V,P)
  [m,n]=size(pi);
  de = Er - pi';
  md = de' * inv(V) * de;
  [mdP,mdQ] = cdfchi("PQ", md, n);
  f=pdfchi(md,n);
  partial=f*(-2*de')*(inv(Omega+P*V*P')*P)';
endfunction


// computeKLIC
// Computes the KL divergence of the posterior
// Inputs
//   pr   - Estimated mean of the posterior distribution
//   pi   - Estimated mean of the prior distribution 
//   V0   - Covariance of the estimated mean of the prior distribution
//   V1   - Covariance of the estimated mean of the posterior distribution
// Outputs
//   kld  - Kullack-Leibler divergence
//   t1   - First term in calculation
//   t2   - second term in calculation
//   t3   - third term in calculation
//   t4   - third term calculated with prior covariance matrix
//
function [kld,t1,t2,t3,t4]=computeKLIC(pr,pi,V0,V1)
  [m,n]=size(pi);
  de = pr - pi;
  t1=log(detr(V1)/detr(V0));
  t2=trace(inv(V1)*V0);
  t3=(de'*inv(V1)*de);
  t4=(de'*inv(V0)*de);
  kld = 0.5*(t1+t2+t3-n);
endfunction

// computeTheil
// Computes compatibility metric for the views vs the prior
// Inputs
//   Omega - Covariance matrix of the views
//   P     - Pick matrix for the views
//   V     - Covariance matrix of prior estimate of mean
//   pi    - Estimated mean of the prior distribution 
//   Q     - Estimated mean of the views
// Outputs
//   meas  - test statistic for Theil's compatibility metric
//   f     - chi square confidence value
//
function [meas,f,partial]=computeTheil(Omega,P,V,pi,Q)
  [m,n]=size(pi);
  denom=inv(Omega + P*V*P')
  m1=P*pi'-Q;
  meas=m1'*denom*m1;
  f=pdfchi(meas,n);
  partial=f*(-2*inv(Omega + P*V*P')*(P*pi'-Q));
endfunction


// Main driver program to run all the scenarios from the paper
// Prints output to the Scilab console that should roughly match
// the tables in paper.

// Display mode
mode(0);
clc;

// Display warning for floating point exception
ieee(1);

// Take the values from He & Litterman, 1999.
weq = [0.016,0.022,0.052,0.055,0.116,0.124,0.615];
C = [ 1.000 0.488 0.478 0.515 0.439 0.512 0.491;
      0.488 1.000 0.664 0.655 0.310 0.608 0.779;
      0.478 0.664 1.000 0.861 0.355 0.783 0.668;
      0.515 0.655 0.861 1.000 0.354 0.777 0.653;
      0.439 0.310 0.355 0.354 1.000 0.405 0.306;
      0.512 0.608 0.783 0.777 0.405 1.000 0.652;
      0.491 0.779 0.668 0.653 0.306 0.652 1.000];
Sigma = [0.160 0.203 0.248 0.271 0.210 0.200 0.187];
refPi = [0.039 0.069 0.084 0.090 0.043 0.068 0.076];
assets=['Australia';'Canada   ';'France    ';'Germany   ';'Japan    ';'UK       ';'USA      '];
labels=['q        ';'omega/tau';'lambda   '];
V = (Sigma' * Sigma) .* C;
 
// Risk tolerance of the market from the paper (page 10)
delta= 2.5;

// Coefficient of uncertainty in the prior estimate of the mean
// From footnote (8) on page 11
tau = 0.05;

// Define view 1
// Germany will outperform the other European markets by 5%
// Market cap weight the P matrix
// Results should match Table 4, Page 21
P1 = [0 0 -.295 1.00 0 -.705 0 ];
Q1 = [0.05];
P=P1;
Q=Q1;
// Define view 2
// Canadian Equities will outperform US Equities by 3%
// Results should match Table 5, Page 22
P2 = [0.0 1.0 0.0 0.0 0.0 0.0 -1.0 ];
Q2 = [0.04];
P = [P1; P2];
Q = [Q1; Q2];
Omega = P * tau * V * P' .* eye(2,2);
[Er,pi,w,postV,lambda] = blacklitterman(delta, weq, V, tau, P, Q, Omega);
t=[100*Q'; diag(Omega)'/tau; lambda'];
output=hypermat([size(assets,1),5]);
format('v',5);
output=[assets string(100*P') string(100*Er) string(100*w)];
bottom=hypermat([3,3]);
format('v',6);
bottom=[labels string(t)];
printf("\nView 1 & 2\n");
mprintf("Country\t\t    P\t\t mu\t w*\n");
mprintf("%s\t%s\t%s\t%s\t%s\n", output);
mprintf("%s\t%s\t%s\n",bottom);
printf("\n");

// Now compute the consistency measure from Fusai and Meucci
[md,mdP,partial]=computeCons(pi,Er,tau*V,P);
printf("\nF&M Consistency and Sensitivities\n %g %g %g %g\n", md, mdP, partial(1), partial(2));

// Compute the TEV
dw=(w-weq');
tev=sqrt(dw'*V*dw);
dq=(tau/delta)*P'*inv((P*tau*V*P') + Omega);
ptev=((V*dw)/tev)'*dq;
printf("TEV and Sensitivities: %g %g %g\n", tev, ptev(1), ptev(2));

// Compute the Kullback-Leibler Information Criteric
klic = - w'*log(w./weq');
mprintf("KLIC (weights) = %g\n", klic);

// Try klic on the returns
klic = - Er' * log(Er./pi');
mprintf("KLIC (returns)  = %g\n", klic);

// KL Divergence between two multi-variate normals
V1=postV;
V0=tau*V;
[kld,t1,t2,t3,t4]=computeKLIC(Er,pi',V0,V1);
mprintf("KLD = %g\n", kld);

[meas,f,partial]=computeTheil(Omega,P,tau*V,pi,Q);
mprintf("Theils measure: %g, T-test %g, Sens: %g %g\n",meas,f,partial(1),partial(2));

[w'*Er,sqrt(w'*(V+V1)*w),w'*Er/sqrt(w'*(V+V1)*w)]

// Raise our confidence in the views by reducing their variance
Omega=Omega*0.25;
[Er,pi,w,postV,lambda] = blacklitterman(delta, weq, V, tau, P, Q, Omega);
t=[100*Q'; diag(Omega)'/tau; lambda'];
output=hypermat([size(assets,1),5]);
format('v',5);
output=[assets string(100*P') string(100*Er) string(100*w)];
bottom=hypermat([3,3]);
format('v',6);
bottom=[labels string(t)];
printf("\nView 1 & 2, View 2 more bullish\n");
mprintf("Country\t\t    P\t\t mu\t w*\n");
mprintf("%s\t%s\t%s\t%s\t%s\n", output);
mprintf("%s\t%s\t%s\n",bottom);
printf("\n");

// Now compute the consistency measure from Fusai and Meucci
[md,mdP,partial]=computeCons(pi,Er,tau*V,P);
printf("\nMahalanobis Distance and View Sensitivities\n %g %g %g %g\n", md, mdP, partial(1), partial(2));

// Compute the TEV
dw=(w-weq');
tev=sqrt(dw'*V*dw);
dq=(tau/delta)*P'*inv((P*tau*V*P') + Omega);
ptev=((V*dw)/tev)'*dq;
printf("TEV and sensitivities: %g %g %g\n", tev, ptev(1), ptev(2));

// Compute the Kullback-Leibler Information Criteria
klic = - w'*log(w./weq');
mprintf("KLIC (weights) = %g\n", klic);

// Try klic on the returns
klic = - Er' * log(Er./pi');
mprintf("KLIC (returns) = %g\n", klic);

// KL Divergence between two multi-variate normals
V1=postV;
V0=tau*V;
[kld,t1,t2,t3,t4]=computeKLIC(Er,pi',V0,V1);
mprintf("KLD = %g\n", kld);

[meas,f,partial]=computeTheil(Omega,P,tau*V,pi,Q);
mprintf("Theils measure: %g, T-test %g, Sens: %g %g\n",meas,f,partial(1),partial(2));

[w'*Er,sqrt(w'*(V+V1)*w),w'*Er/sqrt(w'*(V+V1)*w)]

// Lower our confidence by raising the variance
Omega=Omega*16;
[Er,pi,w,postV,lambda] = blacklitterman(delta, weq, V, tau, P, Q, Omega);
t=[100*Q'; diag(Omega)'/tau; lambda'];
output=hypermat([size(assets,1),5]);
format('v',5);
output=[assets string(100*P') string(100*Er) string(100*w)];
bottom=hypermat([3,3]);
format('v',6);
bottom=[labels string(t)];
printf("\nView 1 & 2, View 2 more bullish\n");
mprintf("Country\t\t    P\t\t mu\t w*\n");
mprintf("%s\t%s\t%s\t%s\t%s\n", output);
mprintf("%s\t%s\t%s\n",bottom);
printf("\n");

// Now compute the consistency measure from Fusai and Meucci
[md,mdP,partial]=computeCons(pi,Er,tau*V,P);
printf("\nMahalanobis Distance and View Sensitivities\n %g %g %g %g\n", md, mdP, partial(1), partial(2));

// Compute the TEV
dw=(w-weq');
tev=sqrt(dw'*V*dw);
dq=(tau/delta)*P'*inv((P*tau*V*P') + Omega);
ptev=((V*dw)/tev)'*dq;
printf("TEV and sensitivities: %g %g %g\n", tev, ptev(1), ptev(2));

// Compute the Kullback-Leibler Information Criteria
klic = - w'*log(w./weq');
mprintf("KLIC (weights) = %g\n", klic);

// Try klic on the returns
klic = - Er' * log(Er./pi');
mprintf("KLIC (returns) = %g\n", klic);

// KL Divergence between two multi-variate normals
V1=postV;
V0=tau*V;
[kld,t1,t2,t3,t4]=computeKLIC(Er,pi',V0,V1);
mprintf("KLD = %g\n", kld);

[meas,f,partial]=computeTheil(Omega,P,tau*V,pi,Q);
mprintf("Theils measure: %g, T-test %g, Sens: %g %g\n",meas,f,partial(1),partial(2));

;[w'*Er,sqrt(w'*(V+V1)*w),w'*Er/sqrt(w'*(V+V1)*w)]

//
// End of bl_metrics.sci