To program a simulation in Matlab that is meant to assist in retirement planning.  Make the following assumptions:

a. You have 30 years until retirement. 

b. Your income is currently $100,000 per year.

c.   It will increase at the rate of inflation (r_inf) which is 3% per year (annualized) (so next year, you’ll earn $103,000).

            d. During each of the next 30 years, you contribute x% of this income to a retirement account.  For simplicity, assume that your contribution is made in a lump sum at the beginning of each year.

            e.   You begin with an asset allocation of 60% stocks and 40% (long-term) bonds.  During each subsequent year, you reduce the asset allocation to stocks, so that after T years, you have 60-T % in stocks and 40+T % in bonds.  (You’ll rebalance your portfolio at the beginning of each year by selling stocks or bonds to ensure that you have the appropriate allocation for that year.)

            f. During each of the 30 years after retirement (years 31-60), you withdraw 70% of your current income, adjusted for inflation.  (For instance, in year 40, you withdraw .7*100,000*(1.03)^40.)

            g. Assume that on average, stocks return 10% per year, while bonds return 6.5%.  Furthermore, the standard deviation of the log-return of stocks over one year is 20%, while the standard deviation of the log-return of bonds is 8%.  The correlation between stock and bond log-returns is 0.2.  (Through time, this correlation is actually very unstable, but more often than not, it is positive.)

The simulation will entail 60 annual timesteps.   What percentage of your income (x) must you contribute to your retirement account to ensure that there is less than a 5% chance that you’ll run out of money during your 30 years of retirement?  What percentage must you contribute to ensure that the probability of exhausting your savings is less than 1%

ms=.1;

ss=.2;

mb=.065;

sb=.08;

rho=.2;

NumTimeSteps=60;

NumPaths=50000;

Initws=.6;

InitAssetValue=100000;

InflaRate=.03;

TimetoRetirement=30;

x=.5;

function[rs,rb]=AssetReturnsSimulation(ms,ss,mb,sb,rho,NumTimeSteps,NumPaths)

dt=1;

y1=randn(NumTimeSteps,NumPaths);

x2=randn(NumTimeSteps,NumPaths);

y2=rho*y1+sqrt(1-rho^2)*x2;

MuStock=log(1+ms)-ss^2/2;

MuBond=log(1+mb)-sb^2/2;

MS=repmat(MuStock,NumTimeSteps,NumPaths);

MB=repmat(MuBond,NumTimeSteps,NumPaths);

rs=exp(MS+ss*sqrt(dt)*y1);

rb=exp(MB+sb*sqrt(dt)*y2);

end 

function [PortfolioReturns]=PortfolioReturnsCalculation(Initws,rs,rb,NumTimeSteps,NumPaths)

ws=zeros(NumTimeSteps,1);

wb=zeros(NumTimeSteps,1);

for t=1:NumTimeSteps

    ws(t)=Initws-t/100;

    wb(t)=1-ws(t);

end

Ws=repmat(ws,1,NumPaths);

Wb=repmat(wb,1,NumPaths);

PortfolioReturns=Ws.*rs+Wb.*rb;

end

function [RetPlan]=CashInflowPattern(InitAssetValue,InflaRate,x,PortfolioReturns,TimetoRetirement,NumPaths)

income=zeros(TimetoRetirement,1);

for t=1:30

    income(t)=InitAssetValue*(1+InflaRate)^t;

end

PR1=PortfolioReturns(1:TimetoRetirement,:);

Income=repmat(income,TimetoRetirement,NumPaths);

CashIn=zeros(TimetoRetirement,NumPaths);

CashIn(1,:)=Income(1,:).*PR1(1,:);

for k=2:30

    CashIn(k,:)=(CashIn(k-1,:)+Income(k,:)).*PR1(k,:); 

end

X=repmat(x,TimetoRetirement,NumPaths);

RetPlan=X.*CashIn;

end

function [RetAcct]=CashOuflowPattern(NumTimeSteps,NumPaths,InitAssetValue,InflaRate,PortfolioReturns,TimetoRetirement,RetPlan)

NumWithdraw=NumTimeSteps-TimetoRetirement;

PR2=PortfolioReturns(TimetoRetirement+1:NumTimeSteps,:);

withdrawal=zeros(NumTimeSteps,1);

for t=TimetoRetirement+1:NumTimeSteps

    withdrawal(t)=0.7*InitAssetValue*(1+InflaRate)^(t);

end

Withdrawal=repmat(withdrawal(TimetoRetirement+1:NumTimeSteps,:),1,NumPaths);

RetAcct(1,:)=(RetPlan(TimetoRetirement,:)-Withdrawal(1,:)).*PR2(1,:);

for t=2:TimetoRetirement;

    RetAcct(t,:)=(RetAcct(t-1,:)-Withdrawal(t,:)).*PR2(t,:);

end

function [RunOutofMoneyPercentage]=Evolve(ms,ss,mb,sb,rho,Initws,InitAssetValue,InflaRate,x,TimetoRetirement,NumTimeSteps,NumPaths)

[rs,rb]=AssetReturnsSimulation(ms,ss,mb,sb,rho,NumTimeSteps,NumPaths);

[PortfolioReturns]=PortfolioReturnsCalculation(Initws,rs,rb,NumTimeSteps,NumPaths);

[RetPlan]=CashInflowPattern(InitAssetValue,InflaRate,x,PortfolioReturns,TimetoRetirement,NumPaths);

[RetAcct]=CashOuflowPattern(NumTimeSteps,NumPaths,InitAssetValue,InflaRate,PortfolioReturns,TimetoRetirement,RetPlan);

Count=zeros(30,1);

for T=1:30;

    Count(T,1)=length(find(RetAcct(T,:)<0));

end

RunOutofMoneyPercentage=max(Count)/length(RetAcct(T,:));

end

function[x]=BisectionAlgorithm(ms,ss,mb,sb,rho,Initws,InitAssetValue,InflaRate,TimetoRetirement,NumTimeSteps,NumPaths,Percentage,Precision)

high=1;

low=0;

x=(high+low)/2;

RunOutofMoneyPercentage=Evolve(ms,ss,mb,sb,rho,Initws,InitAssetValue,InflaRate,x,TimetoRetirement,NumTimeSteps,NumPaths);

while abs(RunOutofMoneyPercentage-Percentage)>Precision

if RunOutofMoneyPercentage>Percentage

        low=x;

else        high=x;

end

x=(high+low)/2;

RunOutofMoneyPercentage=Evolve(ms,ss,mb,sb,rho,Initws,InitAssetValue,InflaRate,x,TimetoRetirement,NumTimeSteps,NumPaths);

end

end

Command Window:

Percentage=.05;

Precision=.00002;

[x]=BisectionAlgorithm(ms,ss,mb,sb,rho,Initws,Initwb,InitAssetValue,InflaRate,TimetoRetirement,NumTimeSteps,NumPaths,Percentage,Precision)

x =0.4293

Percentage=.01;

[x]=BisectionAlgorithm(ms,ss,mb,sb,rho,Initws,Initwb,InitAssetValue,InflaRate,TimetoRetirement,NumTimeSteps,NumPaths,Percentage,Precision)

x =0.5677