// Black-Litterman example code for SciLab (method_idzorek.sci)
// Copyright (c) Jay Walters, blacklitterman.org, 2008.
//
// 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
// "A STEP-BY-STEP GUIDE TO THE BLACK-LITTERMAN MODEL, Incorporating
// user-specified confidence levels" by Thomas Idzorek. You can find
// a copy of this paper at the following url.
//   http://faculty.fuqua.duke.edu/~charvey/Teaching/BA453_2006/Idzorek_onBL.pdf
//
// 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
//

// bl_omega
//   This function computes the Black-Litterman parameters Omega from
//   an Idzorek confidence.
// Inputs
//   conf   - Idzorek confidence specified as a decimal (50% as 0.50)
//   P      - Pick matrix for the view
//   Sigma  - Prior covariance matrix
// Outputs
//   omega  - Black-Litterman uncertainty/confidence parameter
//
function [omega] = bl_omega(conf, P, Sigma)
  alpha = (1 - conf) / conf;
  omega = alpha * P * Sigma * P';
endfunction

// bl_impliedFromW
//   This function computes the implied confidence the tilt of each
//   asset from the equilibrium.
//  Inputs
//    w     - Posterior weight (tilted)
//    weq   - Equilibrium weight
//    w100  - Weight under 100% confidence in the views
//
function [conf] = bl_impliedFromW(w, weq, w100)
  conf = (w - weq)./(w100 - weq);
  [m,n]=size(w);
  conf = conf(1:m,1);
endfunction

// bl_implied
//   This function computes the implied value of the Idzorek confidence
//   given the Black-Litterman parameter Omega.
//
function [conf] = bl_implied(Omega, P, Sigma)
  t = P * Sigma * P';
  alpha = Omega / t;
  conf = 1 / (1 + alpha);
endfunction

// bl_w100
//   This function computes the weights with 100% certainty
// Inputs
//   delta  - Risk tolerance from the equilibrium portfolio
//   weq    - Weights of the assets in the equilibrium portfolio
//   sigma  - Prior covariance matrix
//   P      - Pick matrix for the view(s)
//   Q      - Vector of view returns
// Outputs
//
function [w, er] = bl_w100(delta, weq, sigma, P, Q)
  // Reverse optimize and back out the equilibrium returns
  // This is formula (12) page 6.
  pi = weq * sigma * delta;
  // Compute posterior estimate of the mean
  // This is a simplified version of formula (8) on page 4.
  er = pi' + sigma * P' * inv(P * sigma * P') * (Q - P * pi');
  // Compute posterior weights based on uncertainty in mean
  w = (er' * inv(delta * sigma))';
endfunction

// bl_idzorek
//   This function performs the Black-Litterman blending of the prior
//   and the views into a new posterior estimate of the returns using
//   Idzorek's method.
// Inputs
//   delta  - Risk tolerance from the equilibrium portfolio
//   weq    - Weights of the assets in the equilibrium portfolio
//   sigma  - Prior covariance matrix
//   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
//   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, w, lambda] = bl_idzorek(delta, weq, sigma, P, Q, Omega)
  // Reverse optimize and back out the equilibrium returns
  // This is formula (12) page 6.
  pi = weq * sigma * delta;
  // Compute posterior estimate of the mean
  // This is a simplified version of formula (8) on page 4.
  er = pi' + sigma * P' * inv(P * sigma * P' + Omega) * (Q - P * pi');
  // Compute posterior weights based on uncertainty in mean
  w = (er' * inv(delta * sigma))';
  // 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' - weq)';
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 Idzorek, 2004.
weq = [ .193400 .261300 .120900 .120900 .013400 .013400 .241800 .034900 ];
Sigma = [ .001005 .001328 -.000579 -.000675 .000121 .000128 -.000445 -.000437  ;
 .001328 .007277 -.001307 -.000610 -.002237 -.000989 .001442 -.001535  ;
 -.000579 -.001307 .059852 .027588 .063497 .023036 .032967 .048039  ;
 -.000675 -.000610 .027588 .029609 .026572 .021465 .020697 .029854  ;
 .000121 -.002237 .063497 .026572 .102488 .042744 .039443 .065994  ;
 .000128 -.000989 .023036 .021465 .042744 .032056 .019881 .032235  ;
 -.000445 .001442 .032967 .020697 .039443 .019881 .028355 .035064  ;
 -.000437 -.001535 .048039 .029854 .065994 .032235 .035064 .079958 ];
refPi = [0.0008 0.0067 0.0641 0.0408 0.0743 0.0370 0.0480 0.0660 ];
assets=['US Bonds  ';'Intl Bonds';'US Lg Grth';'US Lg Value';'US Sm Grth';
        'US Sm Value';'Intl Dev Eq';'Intl Emg Eq'];
labels=['q        ';'omega/tau';'lambda   '];
 
// Risk tolerance of the market from the paper (note on page 4)
delta= 3.07;

// Define view 1
// International Developed Equity will have an excess return of 5.25%
// with a confidence of 25%.
P1 = [0 0 0 0 0 0 1 0];
Q1 = [0.0525];

// Define view 2
// International Bonds will outperform US Bonds by 0.0025 with a
// confidence of 50%.
P2 = [-1 1 0 0 0 0 0 0];
Q2 = [0.0025];

// Define View 3
// US Large and Small Growth will outperform US Large and Small Value
// by 0.02 with a confidence of 65%.
P3 = [0 0 0.90 -0.90 0.10 -0.10 0 0];
Q3 = [0.02];

// The pick or focus matrix
P=[P1;P2;P3];
// Make a filter matrix, it is the identity matrix for any
// asset in any view, 0 otherwise
[n,m]=size(P);
F=hypermat([m,m]);
for i = 1:m,
  if (nnz(P(1:n,i))>0) then
    F(i,i)=1;
  end,
end;
Q=[Q1;Q2;Q3];
Omega=diag(diag(P*Sigma*P'));
[Er,w,lambda] = bl_idzorek(delta, weq, Sigma, P, Q, Omega);
t=[100*Q'; diag(Omega)'; lambda'];
output=hypermat([size(assets,1),4]);
format('v',5);
output=[assets string(100*P') string(100*Er) string(100*w)];
bottom=hypermat([3,2]);
format('v',6);
bottom=[labels string(t)];
printf("\nAll Views Default (50%%) Confidence\n");
mprintf("Asset Class\t\tP\t mu\t w*\n");
mprintf("%s\t%s\t%s\t%s\t%s\t%s\n", output);
mprintf("%s\t%s\t%s\t%s\n",bottom);
printf("\n");

// So compute the combined confidence levels, each view has
// also has confidence levels, but they are all 50% to start
[w100,e100]=bl_w100(delta, weq, Sigma, P, Q);
conf=bl_impliedFromW(w, weq', w100);
printf("Aggregate Confidence in asset tilts\n");
(F*conf)'

// The separate implied confidence levels
for i = 1:3,
  conf(i) = bl_implied(Omega(i,i),P(i,1:m),Sigma);
end;
printf("Implied confidence in individual views\n");
conf(1:3)'

// Now get the right confidence levels and put them into Omega
Omega = [bl_omega(0.25, P1, Sigma) 0 0;
        0 bl_omega(0.50, P2, Sigma) 0;
        0 0 bl_omega(0.65, P3, Sigma)];
[Er,w,lambda] = bl_idzorek(delta, weq, Sigma, P, Q, Omega.*eye(Omega));
t=[100*Q'; diag(Omega)'; lambda'];
output=hypermat([size(assets,1),4]);
format('v',5);
output=[assets string(100*P') string(100*Er) string(100*w)];
bottom=hypermat([3,2]);
format('v',6);
bottom=[labels string(t)];
printf("\nAll Views Specified Confidences\n");
mprintf("Asset Class\t\tP\t mu\t w*\n");
mprintf("%s\t%s\t%s\t%s\t%s\t%s\n", output);
mprintf("%s\t%s\t%s\t%s\n",bottom);
printf("\n");

// Compute the aggregated confidence levels
conf=bl_impliedFromW(w, weq', w100);
printf("Aggregate Confidence in asset tilts\n");
(F*conf)'

// The separate implied confidence levels
for i = 1:3,
  conf(i) = bl_implied(Omega(i,i),P(i,1:m),Sigma);
end;
printf("Implied confidence in individual views\n");
conf(1:3)'

//
// End of bl_idzorek.sci