Code covered by the BSD License  

Highlights from
Modeling Variable Annuities with MATLAB

image thumbnail

Modeling Variable Annuities with MATLAB

by

 

This demo shows how to price variable annuity product (Guaranteed Minimum Withdrawal Benefit)

GMWB
classdef GMWB
    
    properties
        riskFreeRate = 0.05; % Default ANNUAL interest rate
        annualWithdrawalRate = 0; % Default annual withrawal Rate
        annualFee = 0; % Default annual fee
        EquityPortfolio = [];
        BondPortfolio = [];
        nPeriods;
        assetPrice = [];
    end
        properties (Access = private)
        nTrials = 1000; % Default no. of trials for Monte Carlo simulations
        EquitySAVal;
        EquityPayoutVal;
        FeeVal;
    end
    
    properties(Dependent = true)
        rfr;   % rfr for stochastic rate
        
    end
    
    methods
        
        function this = GMWB(annualWithdrawalRate, annualFee, rfRate, EP, BP, Periods)
            if ~isempty(annualWithdrawalRate), this.annualWithdrawalRate = annualWithdrawalRate; end
            if ~isempty(annualFee), this.annualFee = annualFee; end
            if ~isempty(rfRate), this.riskFreeRate = rfRate; end
            if ~isempty(EP)
                this.EquityPortfolio = getEquityData(EP.Ticker, EP.FromDate, EP.ToDate, EP.Period);
                this.EquityPortfolio.Holdings = EP.Holdings;
            end
            if ~isempty(BP)
                % To be implemented later
                this.BondPortfolio = BP;
            end
            if ~isempty(Periods)
                this.nPeriods = Periods;
            else
                this.nPeriods = ceil(1/this.annualWithdrawalRate);
            end
            this.assetPrice = this.simAssetPrice();
        end
        
        function [cost, fee, probRuin] = price(this, varargin)
            [this.EquitySAVal, this.EquityPayoutVal, this.FeeVal] = this.calcEquitySAPayoutFeeVal();
            cost = pvvar(mean(this.EquityPayoutVal, 2), this.rfr);
            fee = pvvar(mean(this.FeeVal, 2), this.rfr);
            probRuin = this.calcProbRuin();
      
            if size(varargin,2) == 1
                switch lower(varargin{1})
                    case 'display'
                        % Print figure to Windows Clipboard
                        if isdeployed
                            figure('Visible', 'off');
                        else
                            disp(['Initial Contract Value = ', num2str(this.EquitySAVal(1, 1)), ...
                                '   Insurance cost = ' num2str(cost), ...
                                '   Fee = ' num2str(fee), ...
                                '   Percentage = ' num2str(cost/this.EquitySAVal(1, 1)*100), ...
                                '   Prob. of Ruin = ' num2str(probRuin)]);
                            figure;
                        end
                        
                        set(gcf,'Position',[100 100 800 500]);
                        hold on;
                        [AX, H1, H2] = plotyy(1:this.nPeriods+1, mean(this.EquitySAVal, 2), 1:this.nPeriods+1, [mean(this.EquityPayoutVal, 2) mean(this.FeeVal, 2)], 'plot', 'bar');
                        set(AX(1), 'FontSize', 16);
                        set(AX(1), 'FontName', 'Calibri');
                        set(AX(1), 'XLim', [1 this.nPeriods+2]);
                        set(get(AX(1),'Ylabel'),'String','Sub-Account Value');
                        
                        set(AX(2), 'FontSize', 16);
                        set(AX(2), 'FontName', 'Calibri');
                        set(AX(2), 'XLim', [1 this.nPeriods+2]);
                        set(get(AX(2),'Ylabel'),'String','Cost Incurred & Fee Collected');
                        
                        xlabel('Year to Expiry');
                        title('Simulated GMWB Results');
                        %legend([H1 H2], 'Sub-Account Value', 'Payout', 'Fee Collected', 'Location', 'North');
                        grid on;
                        hold off;
                        
                        if isdeployed
                            print -dmeta
                            close(gcf)
                        end
                    otherwise
                        % do nothing
                end
            end
        end
        
        function val = simAssetPrice(this)
            [expReturn, sigma, correlation] = this.calcRetSigCorr();
            StartAssetPrice = this.EquityPortfolio.TimeSeries(end, :)';
            
            
            %             %%%%%%%%%% For Test
            %             expReturn = 0.06;
            %             sigma = 0.18;
            %             correlation = 1;
            %             %%%%%%%%%% End Test
            
            
            obj = gbm(diag(expReturn), diag(sigma), 'Correlation', correlation, 'StartState', StartAssetPrice);
            
            dt       =   1;
            %strm = RandStream('mt19937ar','Seed',22814);
            strm = RandStream('mt19937ar');
            RandStream.setDefaultStream(strm);
            [val, ~] = obj.simulate(this.nPeriods, 'DeltaTime', dt, 'nTrials', this.nTrials);
            
            val(val < 0) = 0;
        end
        
        function [SAVal, PayoutVal, FeeVal] = calcEquitySAPayoutFeeVal(this)
            % This function post-process simulated asset price
            SAVal = nan(this.nPeriods+1, this.nTrials);
            PayoutVal = nan(this.nPeriods+1, this.nTrials);
            FeeVal = nan(this.nPeriods+1, this.nTrials);
            
            for simIdx = 1:this.nTrials
                curHoldings = this.EquityPortfolio.Holdings;
                curPrice = this.assetPrice(1, 1:end, simIdx);
                SAVal(1, simIdx) =  curPrice * curHoldings;
                PayoutVal(1, simIdx) = 0;
                FeeVal(1, simIdx) = 0;
                guaranteedWithdrawal = curPrice * curHoldings * this.annualWithdrawalRate;
                
                %                 %%%%%%%%% begin test
                %                 guaranteedWithdrawal = 10;
                %                 %%%%%%%%% end test
                
                for i = 2:this.nPeriods+1
                    curPrice = this.assetPrice(i, 1:end, simIdx);
                    curPortVal = curPrice * curHoldings;
                    totalWithdrawal = curPortVal * this.annualFee + guaranteedWithdrawal;
                    FeeVal(i, simIdx) = curPortVal * this.annualFee;
                    
                    [curHoldings, shortfall] = this.updateHoldings(curPrice, curHoldings, totalWithdrawal, 'largest');
                    
                    PayoutVal(i, simIdx) = shortfall;
                    SAVal(i, simIdx) = curPrice * curHoldings;
                end
            end
        end
        
        function [SAVal, PayoutVal] = calcBondSAPayoutVal(this)   %#ok
            % To be implemented
            SAVal = 0;
            PayoutVal = 0;
        end
        
        function [expReturn, sigma, correlation] = calcRetSigCorr(this)
            returns     = price2ret(this.EquityPortfolio.TimeSeries);
            expReturn   = mean(returns);
            sigma       = std(returns);
            correlation = corrcoef(returns);
            
            % Annualize monthey, weekly and daily returns and correlations
            if this.EquityPortfolio.Period == 'm'
                [expReturn, expCov] = arith2geom(mean(returns), cov(returns), 12);
                sigma = diag(sqrt(expCov));
                correlation = corrcov(expCov);
                expReturn = expReturn';
            elseif this.EquityPortfolio.Period == 'w'
                % Not implemented
            elseif this.EquityPortfolio.Period == 'd'
                % Not implemented
            else
                % Do nothing
            end
        end
        
        function val = get.rfr(this)
            % Change this function later to implement stochastic rate
            val = this.riskFreeRate;
        end
    end
    
    methods (Access = private)
        function EP = getEquityData(Ticker, FromDate, ToDate, Period)
            conn = yahoo;
            if isconnection(conn)
                for i = 1:length(Ticker)
                    d = fetch(conn, Ticker(i), 'Close', FromDate, ToDate, Period);
                    sData(:, i) = d(:, 2);                                            %#ok
                end
                EP.TimeSeries = flipud(sData);
                EP.Ticker = Ticker;
                EP.FromDate = FromDate;
                EP.ToDate = ToDate;
                EP.Period = Period;
            else
                d = load('stock.mat');
                EP = d.EP;
            end
        end
        
        function [holdings, shortfall] = updateHoldings(this, curPrice, curHoldings, totalWithdrawal, method)
            shortfall = 0;
            holdings = curHoldings;
            curPrice = curPrice';
            curPortVal = curHoldings .* curPrice;
            
            if any(curPortVal > 0)
                curPortVal(curPortVal > 0) = curPortVal(curPortVal > 0) - totalWithdrawal/sum(curPortVal > 0);
                newPortVal = curPortVal;
                if sum(newPortVal) > 0      % Still have fund left in the account
                    if all(newPortVal > 0)  % All stocks have positive values
                        shortfall = 0;
                        holdings = newPortVal ./ curPrice;
                        holdings(~isfinite(holdings)) = 0;
                    else   % Some stocks are negative, meaning we have lost more than we invested
                        additionalLoss = -sum(newPortVal(newPortVal < 0));
                        pVal = newPortVal;
                        pVal(pVal < 0) = 0;
                        while additionalLoss > 0
                            switch lower(method)
                                case 'random'
                                    allIdx = find(pVal > 0);
                                    index = allIdx(randi(numel(allIdx), 1));
                                case 'largest'
                                    [~, index]  = max(pVal);   % Unload from largest holding first
                                case 'smallest'
                                    index = find((pVal > 0) & (pVal == min(pVal(pVal > 0)))); % Unload from smallest holding first
                                otherwise
                                    index = find((pVal > 0) & (pVal == min(pVal(pVal > 0))));
                            end
                            
                            if pVal(index) > additionalLoss
                                pVal(index) = pVal(index) - additionalLoss;
                                additionalLoss = 0;
                            else
                                additionalLoss = additionalLoss - pVal(index);
                                pVal(index) = 0;
                            end
                        end
                        shortfall = 0;
                        holdings = pVal ./ curPrice;
                        holdings(~isfinite(holdings)) = 0;
                    end
                else
                    shortfall = -sum(newPortVal);
                    holdings(:) = 0;
                end
            else
                shortfall = totalWithdrawal;
                holdings(:) = 0;
            end
        end
        
        function probRuin = calcProbRuin(this)
            probRuin = sum(min(this.EquitySAVal) == 0) / size(this.EquitySAVal, 2);
        end
    end
    
end

Contact us