// global-equilibrium.sci
// This scilab (should be MATLAB compatible as well) script takes the numbers from the
// examples in Litterman (2003) Chapter 6 where he computes the Global Equilibrium
// asset returns and weights given the covariance matrix and a suitably rich set of
// equilibrium constraints.
//
// Jay Walters
// jwalters@blacklitterman.org
// 14 May 2008

// Clear the Scilab command window
clc;
// Set the display size to accomodate our matrices in one page
nl=lines();
lines(40,132);

// First we try the simple 2 country, 2 asset model

// Final results for returns, we will use these to check against 
muu = [0.04128 0.0323 0.00412];
muj = [0.04038 0.0306 0.00588];

// Final results for weights we will use these to check our results against
du = [0.80 0.20 0.10 -0.10];
wu = [0.80 0.20 0.10];
dj = [0.80 0.20 -0.40 0.40];
wj = [0.80 0.20 0.40];

// The expected final solution vector will be 
X=hypermat(1,[14]);
X(1:3)=muu;
X(4:7)=du;
X(8:10)=muj;
X(11:14)=dj;

// This is the US investors inverse covariance matrix
// This is slightly more accurate and tries to take into account some rounding error
//iSIGMAu = [ 59.2661 -26.093 -0.89834; -26.093 46.43867 -5.546 ; -.89834 -5.546 101.0242 ];
iSIGMAu = [ 59.27 -26.09 -0.90; -26.09 46.44 -5.55 ; -.90 -5.55 101.02 ];
// This is the US investors covariance matrix
SIGMAu = inv(iSIGMAu);

// This is the JP investors inverse covariance matrix
iSIGMAj = [iSIGMAu(1,1) iSIGMAu(1,2) -iSIGMAu(1,3);iSIGMAu(2,1) iSIGMAu(2,2) -iSIGMAu(2,3); -iSIGMAu(3,1) -iSIGMAu(3,2) iSIGMAu(3,3)];
// These are the values in the book (which should match the stuff above)
iSIGMAj = [59.27 -26.09 0.90 ; -26.09 46.44 5.55 ; 0.90 5.55 101.02 ];
// This is hte JP investors covariance matrix
SIGMAj = inv(iSIGMAj);

// Now we compute the terms from Siegel's paradox
// mu($u) - mu(yu) = sigma($xu)
sigmadxu = SIGMAu(3,1);
// mu($j) - mu(yj) = sigma($xj)
sigmajxu = SIGMAu(3,2);

// We can test the final results to make sure they match what Siegel's paradox predicts
muu(1)-muj(1);
muu(2)-muj(2);

// Now on to demonstrate the problem, we are solving a simple linear system Ax=b, so we need to
// specify A and b, x will be the computed answer.

// Now create the A matrix and the b vector, and initialize them to all 0's
A=hypermat([14,14]);
b=hypermat([14]);

// Copy the inverse covariance matrices into the right place in A
lambda=2.0;
A(1:3,1:3)=(1/lambda)*iSIGMAu;
A(1:2,4:5)=-1*eye(2,2);
A(3,7)=-1;
A(4:6,8:10)=(1/lambda)*iSIGMAj;
A(4:6,11:13)=-1*eye(3,3);

// Given that we have lambda = 2, then we should expect our hedging fraction to be 50%

// Now the constraints due to Siegel's paradox
// the difference in us equity returns = the covariance between currency and us equity
A(7,1)=1;
A(7,8)=-1;
b(7)=SIGMAu(3,1);
// The difference in jp equity returns = the covariance between currency and jp equity
A(8,2)=1;
A(8,9)=-1;
b(8)=SIGMAu(3,2);
// The sum of the usd and jpy returns = the variance of returns
A(9,3)=1;
A(9,10)=1;
b(9)=SIGMAu(3,3);

// Now the constraints on total wealth
// us wealth * us fraction in us equity + jp wealth * jp fraction in us equity = us equity market cap
A(10,4)=80;
A(10,11)=20;
b(10)=80;
// us wealth * us fraction in jp equity + jp wealth * jp fraction in jp equity = jp equity market cap
A(11,5)=80;
A(11,12)=20;
b(11)=20;

// Now the constraints on lending
// us wealth * us dollar lending + jp wealth * jp dollar borrowing = 0
A(12,6)=80;
A(12,13)=20;
b(12)=0;

// us wealth * us yen borrowing + jp wealth * jp yen lending = 0
A(13,7)=80;
A(13,14)=20;
b(13)=0;

// And the constraints on total investment
// We only need one to get us 14x14, the sum of us fractions = 1
A(14,4:7)=1;
b(14)=1;

// Now try and solve the system, linsolve solves AX + b = 0 so we need to negate b as
// we computed the system Ax = b above.
x=linsolve(A,-b');

labels=['US Equity Return (USD)';'JP Equity Return (USD)';'JPY Return (USD)';
        'US Investor US Equity';'US Investor JP Equity';'US Investor USD Lending';'US Investor JPY Lending';
        'US Equity Return (JPY)';'JP Equity Return (JPY)';'USD Return (JPY)';
        'JP Investor US Equity';'JP Investor JP Equity';'JP Investor USD Lending';'JP Investor JPY Lending'];

// Print out the computed solution and the one from the example in the book
format('v',8);    // Not sure why this needs to be 8 not 4
printf("2 Country 2 Asset Global Equilibrium Results\n");
for i = 1:14, printf("%s\tx[%02d] = %f (%f)\n", labels(i),i, x(i), x(i)-X(i)), end

// Check the hedging ratios
printf("\n");
printf("US Investor Hedges %f%%\n", -100.0*x(7)/x(5));
printf("JP Investor Hedges %f%%\n", -100.0*x(13)/x(11));

// Now play with the quick way
I1 = [1 0 0;0 1  0;0 0 1];
nSIGMAu=I1*SIGMAu*I1';
I2 = [1 0 0;0 1 0;0 0 -1];
nSIGMAj=I2*SIGMAu*I2';

H1=[1 0 0; 0 1 0;-1 0 -1;0 -1 1];
nI1=[0 0 1 0];
H2=[1 0 0;0 1 0;-1 0 1;0 -1 -1];
nI2=[0 0 0 1];
d1 = nI1' + H1*wu'
d2 = nI2' + H2*wj'

J = H1 * inv(I1)';

np1I2 = [0 0 1]';
mu2 = I2*muu'+ (I2*SIGMAu*I2')*np1I2;
// Matches muj right on, so given mu1 we can get to mu2 easily

W=[0 0 0.80 0.20]';
s=W(3)*d1+W(4)*d2;
s=[0.80 0.20 0 0]';

// According to the book, j should be as follows
jp=W(3)*(1/lambda)*H1*np1I2+W(4)*(1/lambda)*H2*np1I2;
// However in another place it says it should be (which is correct)
j=W(4)*(1/lambda)*H2*np1I2;
tau=W(3)*(1/lambda)+W(4)*(1/lambda);

// Test eqn 6.26
s2=W+tau*J*iSIGMAu*muu' + j;
mu1=(1/tau)*SIGMAu*pinv(J)*(s-W-j);

// reverse engineer j - just to check what we have
s-W-J*iSIGMAu*muu'*tau

wu2=(1/(tau*lambda))*inv(J'*J)*J'*(s-W-j)
wj2=inv(I2)'*wu'+(1/lambda)*np1I2

// Reset display width/height
lines(nl(2),nl(1));