disp('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

disp('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');

disp('Mean:');

disp(hrznInvMu);

disp('Standard Deviation:');

disp(hrznInvStd);

disp('Correlation:');

disp(hrznInvCor);

disp('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');

disp('Mean:');

disp(hrznInvSimsMu);

disp('Standard Deviation:');

disp(hrznInvSimsStd);

disp('Correlation:');

disp(hrznInvSimsCor);

disp('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'});

disp('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');

disp(investorPortfolioWts);