clear
echo on
format compact
%Solution to Problem 23
%In this problem we are asked to maximize
%the return of an initial investment subject
%to having a 97.5% chance of some minimum
%return after 20 years.  This percent figure
%means that the expected value after 20 years
%must be two standard deviations above the
%minimum objective.
pause

%To solve this problem, we need to write a
%routine which calculates the expected return
%and the 97.5% probability floor of the return
%for a given investment ratio.  We do this in
%the attached function yield.m which should be
%saved in a separate file.  We make the assumption
%that the yields produced by the three investment
%vehicles are independent, and that there is no
%covariance from year to year either.  These are
%probably rather unrealistic assumptions, but it
%does make the calculation simpler!
pause

%Because we will invest the money in one of three
%possible ways, there are only two independent variables
%describing the division.  We let x(1) be the fraction
%put in T-bills, and x(2) be the amount put in stocks.
%The remaining fraction 1-x(1)-x(2) is put in the
%Asian market fund.  Because we have a constrained
%optimization problem, we introduce the slack variable
%x(3) to convert this to an equality constraint.  We
%will use the penalty function approach to determine
%the optimum yield subject to the imposed constraint.
pause

%We use the function funmin.m given below, which in turn
%calls the function yield.m, to determine the optimum
%conditions.  We pass the value of the penalty parameter
%P, the initial investment initval, the number of years
%of investment years, and the goal after these years goal
%to the functions as global parameters.  Thus:
global P initval goal years
initval=50000;
goal=150000;
years=20;
pause

%This is our initial guess:
x=ones(3,1)/3;
%and our initial value of P:
P=100;
plot(x(2),1-x(1)-x(2),'o')
axis([0,1,0,1])
xlabel('Stock fund investment')
ylabel('Asia fund investment')
hold on
plot([0,1],[1,0],'r')
drawnow
xt=x+1;
OPTIONS(1)=0;
OPTIONS(2)=1e-5;
while norm(x-xt)>.00001,
 xt=x;
 x=fmins('funmin',x,OPTIONS);
 x=abs(x);%we force x to be positive
 P=P*2;
 plot(x(2),1-x(1)-x(2),'o')
 drawnow
 echo off
end
echo on
pause

%The optimum values are thus:
x=abs(x);
Tbills=x(1)
stocks=x(2)
Asiafund=1-x(1)-x(2)
[ret,floor]=yield(x)
%And the final value of the penalty function
%parameter was:
P
pause

%Note that nothing is being put in T-bills, (you
%can see this because the optimum is on the red
%line) and the funds are divided about evenly between
%the stock fund and the Asia fund.  We had to start
%out with a fairly large value of P in order to
%avoid getting frozen out at the point where
%everything is in the Asia fund and the minimum
%goal constraint is not satisfied.
pause

%Let's do this again for the case where we have
%a shorter term investment:
years=5;
initval=75000;
goal=100000;
P=100;
x=ones(3,1)/3;
plot(x(2),1-x(1)-x(2),'+w')
xt=x+1;
OPTIONS(1)=0;
OPTIONS(2)=1e-5;
while norm(x-xt)>.00001,
 xt=x;
 x=fmins('funmin',x,OPTIONS);
 x=abs(x); %we force x to be positive here
           %and in the function yield.m
 P=P*2;
 plot(x(2),1-x(1)-x(2),'+g')
 drawnow
 echo off
end
echo on
pause

hold off
%The optimum values are now:
x=abs(x);
Tbills=x(1)
stocks=x(2)
Asiafund=1-x(1)-x(2)
[ret,floor]=yield(x)
%And the final value of the penalty function
%parameter was:
P
pause

%Note that the mix has now changed, so that
%we invest mostly in T-bills.  This is why most
%investment fund managers recommend that investors
%gradually move their funds from growth (high yield /
%high risk) investments to more secure investments as
%the time to utilization (e.g., college tuition or
%retirement) approaches.
echo off


__________________________


function y=funmin(x)
%This function returns the quantity to be
%minimized in the optimization routine.  We
%want to maximize the yield subject to some
%minimum floor.  We make use of a slack variable
%and a penalty function:
global P goal

%we require all of the investments to be positive,
%thus:
x=abs(x);
[ret,floor]=yield(x);
g=(floor-goal)/goal-x(3)^2;  %We normalize the magnitude of g to unity.

y=-ret+P*norm(g)^2;


___________________________


function [ret,floor]=yield(x)
%This function returns the expected multiplication
%of an initial investment amount as a function of
%the distribution of the investments, and also
%returns the 97.5% confidence lower bound, or floor,
%of the return as well.  The number of years invested
%and the initial investment are passed in as global
%parameters.

global initval years
%The formula for calculating the expected return is
%simple.  Each year each account increases by an amount
%proportional to the balance at the beginning of the
%year and the yield.  Thus:
x=abs(x); %We force x to be positive!
bond=x(1)*(1+.07)^years;
stock=x(2)*(1+0.1)^years;
asia=(1-x(1)-x(2))*(1+0.15)^years;
ret=bond+stock+asia;

%The error is a bit harder to calculate.  Because the
%variability is pretty large, let's look at the log of
%the return.  We are calculating the gradient of the
%log of the return with respect to the yield in each
%year, thus:
gradt=x(1)*(1.07)^(years-1)*ones(1,years);
grads=x(2)*(1.1)^(years-1)*ones(1,years);
grada=(1-x(1)-x(2))*(1.15)^(years-1)*ones(1,years);
gradret=[gradt,grads,grada]/ret;

%We also need the matrix of covariance of the yields
%in each year of the funds.  If we assume that they are
%all independent, then it is just a diagonal matrix:

varyield=diag([0*ones(1,years),0.1^2*ones(1,years),0.25^2*ones(1,years)]);

%Note that this would be a more complex matrix if we
%had properly modelled the covariance of the yields!

%The variance in the return is thus:
varlogret=gradret*varyield*gradret';

%and the floor of the initial investment is:
floor=initval*exp(log(ret)-2*varlogret^.5);

%If we didn't do it this way, we would have gotten
%into trouble with the high degree of variability in
%the asia fund.  The worst you can do is lose all of the
%money invested in that fund - not losing more than that
%amount!  Working with the log of the returns assures us
%that the floor will be a positive number for any degree
%of variability.  It is also much more accurate, as the
%log of the return is close to being a normally distributed
%variable, while the return itself is not.