Main Content

Bond Portfolio Optimization

This example shows how to construct an optimal portfolio of 10,20 and 30 year treasuries that will be held for a period of one month. Emphasis is placed on the over-all asset allocation process.

  • Step 1: Load Market Data - Historic daily treasury yields downloaded from FRED are loaded.

  • Step 2: Calculate Market Invariants - Daily changes in yield to maturity are chosen as invariants and assumed to be multivariate normal. Due to missing data for the 30 year bonds, an expectation maximization algorithm is used to estimate the mean and covariance of the invariants. The invariant's statistics are projected to the investment horizon.

  • Step 3: Simulate Invariants at Horizon - Due to the high correlation and inherent structure in the yield curves, a principal component analysis is applied to the invariant statistics. Multivariate normal random draws are done in PCA space. The simulations are transformed back into invariant space using the PCA loadings.

  • Step 4: Calculate Distribution of Returns at Horizon - The simulated monthly changes in the yield curve are used to calculate the yield for the portfolio securities at the horizon. This requires interpolating values off of the simulated yield curves since the portfolio securities will have maturities that are one month less than 10, 20 and 30 years. Profit/Loss for each scenario/security are calculated by pricing the treasuries using the simulated and interpolated yields. Simulated linear returns and their statistics are calculated from the prices.

  • Step 5: Optimize Asset Allocation - Quadratic mean/variance optimization is performed on the treasury returns statistics to calculate optimal portfolio weights for 10 points along the Efficient Frontier. The investor preference is to choose the portfolio that is closest to the mean value of possible Sharpe ratios.

Step 1: Load Market Data

Historic Yield To Maturity Data for Series: DGS6MO, DGS1, DGS2, DGS3, DGS5, DGS7, DGS10, DGS20, DGS30 For Dates: Sep 1, 2000 - Sep 1, 2010 Obtained from: http://research.stlouisfed.org/fred2/categories/115 Note: Data is downloaded using Datafeed Toolbox™ using commands like: >> conn = fred; >> data = fetch(conn,'DGS10','9/1/2000','9/1/2010'); Results have been aggregated and stored in a binary MAT file for convinience

disp('Step 1: Load and Visualize Market Data...');
Step 1: Load and Visualize Market Data...
histData = load('HistoricalYTMData.mat');
% Time to maturity for each series
tsYTMMats = histData.tsYTMMats;
% Dates rates were observed
tsYTMObsDates = histData.tsYTMObsDates;
% Observed Rates
tsYTMRates = histData.tsYTMRates;

% Visualize Yield Curves
miny = min(tsYTMRates(:));
maxy = max(tsYTMRates(:));

figure;
h = plot(tsYTMMats,tsYTMRates,'k-o');
axis([0,32,miny,maxy]);
xlabel('Time to Maturity');
ylabel('Yield');
legend('Historic Yield Curve','location','se');
grid on;
for i = 1:50:length(tsYTMObsDates)
    set(h,'ydata',tsYTMRates(i,:));
    title(datestr(tsYTMObsDates(i)));
    pause(0.1);
end

Step 2: Calculate Market Invariants

For market invariants, use the standard: daily changes in yield to maturity for each series. You can estimate their statistical distribution to be multi-variate normal. IID Analysis on each invariant series produces decent results - more so in the "independent" factor than "identical". A more thorough modeling using more complex distributions and/or time series models is beyond the scope of this project. What will need to be accounted for is the estimation of distribution parameters in the the presence of missing data. The 30 year bonds were discontinued for a period between Feb 2002 and Feb 2006, so there are no yields for this time period.

disp('Step 2: Calculate Market Invariants...');
Step 2: Calculate Market Invariants...
% Invariants are assumed to be daily changes in YTM rates
tsYTMRateDeltas = diff(tsYTMRates);

% About 1/3 of the 30 year rates (column 9) are missing from the original
% data set. Rather than throw out all these observations, an expectation
% Maximization routine is used to estimate the mean and covariance of the
% invariants. Default options (NaN skip for initial estimates, etc.) are used.
[tsInvMu,tsInvCov] = ecmnmle(tsYTMRateDeltas);

% Calculate standard deviations and correlations
[tsInvStd,tsInvCorr] = cov2corr(tsInvCov);

% The investment horizon is 1 month. (21 business days between 9/1/2010
% and 10/1/2010). Since the invariants are summable and the means and
% variances of normal distributions are normal, we can project the
% invariants to the investment horizon as follows
hrznInvMu = 21*tsInvMu';
hrznInvCov = 21*tsInvCov;
[hrznInvStd,hrznInvCor] = cov2corr(hrznInvCov);

% Display results
disp('The market invariants projected to the horizon have the following stats');
The market invariants projected to the horizon have the following stats
disp('Mean:');
Mean:
disp(hrznInvMu);
   1.0e-03 *

   -0.5149   -0.4981   -0.4696   -0.4418   -0.3788   -0.3268   -0.2604   -0.2184   -0.1603
disp('Standard Deviation:');
Standard Deviation:
disp(hrznInvStd);
    0.0023    0.0024    0.0030    0.0032    0.0033    0.0032    0.0030    0.0027    0.0026
disp('Correlation:');
Correlation:
disp(hrznInvCor);
    1.0000    0.8553    0.5952    0.5629    0.4980    0.4467    0.4028    0.3338    0.3088
    0.8553    1.0000    0.8282    0.7901    0.7246    0.6685    0.6175    0.5349    0.4973
    0.5952    0.8282    1.0000    0.9653    0.9114    0.8589    0.8055    0.7102    0.6642
    0.5629    0.7901    0.9653    1.0000    0.9519    0.9106    0.8664    0.7789    0.7361
    0.4980    0.7246    0.9114    0.9519    1.0000    0.9725    0.9438    0.8728    0.8322
    0.4467    0.6685    0.8589    0.9106    0.9725    1.0000    0.9730    0.9218    0.8863
    0.4028    0.6175    0.8055    0.8664    0.9438    0.9730    1.0000    0.9562    0.9267
    0.3338    0.5349    0.7102    0.7789    0.8728    0.9218    0.9562    1.0000    0.9758
    0.3088    0.4973    0.6642    0.7361    0.8322    0.8863    0.9267    0.9758    1.0000

Step 3: Simulate Market Invariants at Horizon

The high correlation is not ideal for simulation of the distribution of invariants at the horizon (and ultimately security prices). A principal component decomposition is used to extract orthogonal invariants. This could also be used for dimension reduction, however since the number of invariants is still relatively small, retain all 9 components for more accurate reconstruction. However, missing values in the market data prevents you from estimating directly off of the time series data. Instead, this can be done directly off of the covariance matrix

disp('Step 3: Simulate Market Invariants at Horizon...');
Step 3: Simulate Market Invariants at Horizon...
% Perform PCA decomposition using invariants' covariance
[pcaFactorCoeff,pcaFactorVar,pcaFactorExp] = pcacov(hrznInvCov);

% Keeping all components of pca decompositon
numFactors = 9;

% Create PCA factor covariance matrix
pcaFactorCov = corr2cov(sqrt(pcaFactorVar),eye(numFactors));

% The number of simulations (random draws)
numSim = 10000;

% Fix random seed for reproducible results
stream = RandStream('mrg32k3a');
RandStream.setGlobalStream(stream);

% Take random draws from multi-variate normal distribution with zero mean
% and diagonal covariance
pcaFactorSims = mvnrnd(zeros(numFactors,1),pcaFactorCov,numSim);

% Transform to horizon invariants and calculate statistics
hrznInvSims = pcaFactorSims*pcaFactorCoeff' + repmat(hrznInvMu,numSim,1);
hrznInvSimsMu = mean(hrznInvSims);
hrznInvSimsCov = cov(hrznInvSims);
[hrznInvSimsStd,hrznInvSimsCor] = cov2corr(hrznInvSimsCov);

% Display results
disp('The simulated invariants have very similar statistics to the original invariants');
The simulated invariants have very similar statistics to the original invariants
disp('Mean:');
Mean:
disp(hrznInvSimsMu);
   1.0e-03 *

   -0.4983   -0.5002   -0.4832   -0.4542   -0.4031   -0.3597   -0.2867   -0.2515   -0.1875
disp('Standard Deviation:');
Standard Deviation:
disp(hrznInvSimsStd);
    0.0023    0.0023    0.0030    0.0031    0.0032    0.0031    0.0029    0.0027    0.0026
disp('Correlation:');
Correlation:
disp(hrznInvSimsCor);
    1.0000    0.8527    0.5827    0.5502    0.4846    0.4327    0.3896    0.3197    0.2961
    0.8527    1.0000    0.8227    0.7840    0.7181    0.6603    0.6097    0.5246    0.4896
    0.5827    0.8227    1.0000    0.9646    0.9100    0.8569    0.8048    0.7074    0.6633
    0.5502    0.7840    0.9646    1.0000    0.9507    0.9085    0.8656    0.7757    0.7344
    0.4846    0.7181    0.9100    0.9507    1.0000    0.9721    0.9428    0.8710    0.8319
    0.4327    0.6603    0.8569    0.9085    0.9721    1.0000    0.9726    0.9211    0.8870
    0.3896    0.6097    0.8048    0.8656    0.9428    0.9726    1.0000    0.9552    0.9264
    0.3197    0.5246    0.7074    0.7757    0.8710    0.9211    0.9552    1.0000    0.9753
    0.2961    0.4896    0.6633    0.7344    0.8319    0.8870    0.9264    0.9753    1.0000

Step 4: Calculate Distribution of Security Returns at Horizon

The portfolio will consist of 10, 20, and 30 year maturity treasuries. For simplicity, assume that these are new issues on the settlement date and are priced at market value inferred from the current yield curve. Profit and Loss distributions are calculated by pricing each security along each simulated yield at the horizon and subtracting the purchase price. The horizon prices require nonstandard time to maturity yields. These are calculated using cubic spline interpolation. Simulated linear returns are their statistics that are calculated from the Profit and Loss scenarios.

disp('Step 4: Calculate Distribution of Security Returns at Horizon...');
Step 4: Calculate Distribution of Security Returns at Horizon...
% Purchase and investment horizon dates
settleDate = '9/1/2010';
hrznDate = '10/1/2010';

% The maturity dates for new issue treasuries purchased on the settle date
treasuryMaturities = {'9/1/2020','9/1/2030','9/1/2040'};

% The observed yields for the securities of interes on the settle date
treasuryYTMAtSettle = tsYTMRates(end,7:9);

% Initialize arrays for later use
treasuryYTMAtHorizonSim = zeros(numSim,3);
treasuryPricesAtSettle = zeros(1,3);
treasuryPricesAtHorizonSim = zeros(numSim,3);

% Use actual/actual day count basis with annualized yields
basis = 8;

% Price treasuries at settle date using known yield to maturity
% Note: For simplicity, we are assuming that none of these securities
% include coupon payments. The hope is that although the prices may not be 
% accurate the overall structure/relationships between value will be
% preserved for the asset allocation process.
for j=1:3
    treasuryPricesAtSettle(j) = bndprice(treasuryYTMAtSettle(j),0,settleDate,...
                                         treasuryMaturities(j),'basis',basis);
end

% To price the treasuries at the horizon, we need to know yield to maturity
% at 9 years 11 months, 19 years 11 months and 29 years 11 months for each
% simulation. We approximate these using cubic spline interpolation

% Transform simulated invariants to YTM at horizon
hrznYTMRatesSims = repmat(tsYTMRates(end,:),numSim,1) + hrznInvSims;

hrznYTMMaturities = {'4/1/2011','10/1/2011','10/1/2012','10/1/2013',...
                     '10/1/2015','10/1/2017','10/1/2020','10/1/2030',...
                     '10/1/2040'};

% Convert dates to numeric serial dates                 
x = datenum(hrznYTMMaturities);
xi = datenum(treasuryMaturities);

% For numerical accuracy, shift x values to start at zero
minDate = min(x);
x = x - minDate;
xi = xi - minDate;

% For each simulation and maturity approximate yield near 10,20,30 year
% nodes. Note that the effects of a spline fit vs. linear fit have a
% significant effect on the resulting ideal allocation. This is due
% to significant under-estimation of yield when using a linear fit
% for points just short of the known nodes
for i=1:numSim
	treasuryYTMAtHorizonSim(i,:) = interp1(x,hrznYTMRatesSims(i,:),xi,'spline');
end

% Visualize 1 simulated yield curve with interpolation
figure;
plot(x,hrznYTMRatesSims(1,:),'k-o',xi,treasuryYTMAtHorizonSim(1,:),'ro');
xlabel('Time (days)');
ylabel('Yield');
legend({'Simulated Yield Curve','Interpolated Yields'},'location','se');
grid on;
title('Zoom to See Spline vs. Linear Interpolants');

% Price treasuries at horizon for each simulated yield to maturity
% Same assumptions are being made here as in above call to bndprice
basis = 8*ones(numSim,1);
for j=1:3
    treasuryPricesAtHorizonSim(:,j) = bndprice(treasuryYTMAtHorizonSim(:,j),0,...
                                               hrznDate,treasuryMaturities(j),'basis',basis);
end

% Calculate distribution of linear returns
treasuryReturns = ( treasuryPricesAtHorizonSim - repmat(treasuryPricesAtSettle,numSim,1) )./repmat(treasuryPricesAtSettle,numSim,1);

% Calculate returns statistics
retsMean = mean(treasuryReturns);
retsCov  = cov(treasuryReturns);
[retsStd,retsCor] = cov2corr(retsCov);

% Visualize results for 30 year treasury
figure;
hist(treasuryReturns,100);
title('Distribution of Returns for 10, 20, 30 Year Treasuries');
grid on;
legend({'10 year','20 year','30 year'});

Step 5: Optimize Asset Allocation

Asset allocation is optimized using quadratic programming. Ten optimal portfolios are calculated and their sharpe ratios are calculated. The optimal portfolio based on investor preference is chosen to be the one that is closest to the mean value of the Sharpe ratio

disp('Step 5: Optimize Asset Allocation...');
Step 5: Optimize Asset Allocation...
% Calculate 10 points along the projection of the efficient frontier into
% std/return space.
[portStd,portRet,portWts] = portopt(retsMean,retsCov,10);

% Visualize
figure;
subplot(2,1,1)
plot(portStd,portRet,'k-o');
xlabel('Portfolio Std');
ylabel('Portfolio Return');
title('Efficient Frontier Projection');
legend('Optimal Portfolios','location','se');
grid on;

subplot(2,1,2)
bar(portWts,'stacked');
xlabel('Portfolio Number');
ylabel('Portfolio Weights');
title('Percentage invested in each treasury');
legend({'10 year','20 year','30 year'});

% Sharpe ratio is calculated using 0 risk-free rate since we are investing
% in treasuries
sharpe = portRet./portStd;

% Investor chooses portfolio based on a near mean Sharpe ratio
sharpeTarget = mean(sharpe);
investorChoice = find( min(abs(sharpe-sharpeTarget)) == abs(sharpe-sharpeTarget));
investorPortfolioWts = portWts(investorChoice,:);

disp('Investor percentage allocation in 10,20,30 year treasuries');
Investor percentage allocation in 10,20,30 year treasuries
disp(investorPortfolioWts);
    0.3989    0.3196    0.2815

See Also

| | | | |

Related Examples

More About