%% Run the following script to generate data for a cost surface
close all; clear all; clc;
annualFee = 0;
rfr = 0.05;
annualWithdrawalRate = 0.10;
EP.Ticker = {'MMM', 'AA', 'AXP', 'T'};
%EP.Ticker = {'MMM'};
EP.FromDate = '01/01/2000';
EP.ToDate = '01/01/2010';
EP.Period = 'm';
EP.Holdings = 10*ones(length(EP.Ticker), 1);
Periods = 10:-1:1;
incrHoldings = 0.1:1:2;
cost = nan(length(incrHoldings), length(Periods));
tic;
for i = 1:length(incrHoldings)
EP.Holdings = incrHoldings(i) * ones(length(EP.Ticker), 1);
for j = 1:length(Periods)
va = GMWB(annualWithdrawalRate, annualFee, rfr, EP, [], Periods(j));
cost(i, j) = va.price();
end
end
toc;
%% Set up grids for surface fitting tool
% PLEASE NOTE that data contained in cost.mat is not the same as those
% generated above. No real data was used to generated data in cost.mat. We
% used a theoretical case of one stock with 10% withdrawal rate
%load cost.mat
[initSAGrid, IGWBGrid] = meshgrid(initSA, IGWB);
initSAGrid = reshape(initSAGrid, numel(initSAGrid), 1); % Y-axis
IGWBGrid = reshape(IGWBGrid, numel(IGWBGrid), 1); % X-axis
costGrid = reshape(cost', numel(cost), 1); % Z-axis
%% Open sftool
% Load init_SA as Y, IGWB_grid as X and cost_grid as Z
% Set Axis Limit from 10 to 0
% Fit with Lowess and Quadratic Polynomial
% Generate M-file
[fittedModel, ~] = createSurfaceFit(IGWBGrid, initSAGrid, costGrid);
%% Plot delta
[~, delta] = differentiate(fittedModel, IGWBGrid, initSAGrid);
[fittedModelDelta, ~] = createSurfaceFit(IGWBGrid, initSAGrid, delta);