# 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 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 ```