* suppress some output to keep list file small
$offlisting
$offsymxref
$offsymlist
$offuelxref
$offuellist
$offupper
option sysout = off;
option solprint = off;

** Inputs data from input file
$include 'portfolio.input'

Parameter Return(I)	Annual return on securities - percentages;
 Return(I) = RF + sum(J,IndBeta(I,J)*Factor(J));

Variables
X(I)	Portion invested in security I
Beta(J)	Factor betas of portfolio
Ret	Return on portfolio
TotRisk	Total betas in portfolio
Positive variable X;

Equations
Fract	Must total to 1
Pbeta(J)	Portfolio beta for factor J
Pret	Portfolio return
BetaCons(J)	Beta of portfolio must be less than DesBeta
RetCons	Return of portfolio must be more than DesRet
Risk	Weighted average of portfolio betas;

Fract..	sum(I,X(I)) =e= 1;
Pbeta(J)..	Beta(J) =e= sum(I,X(I)*IndBeta(I,J));
Pret..	Ret =e= sum(I,X(I)*Return(I));
BetaCons(J).. Beta(J) =l= DesBeta(J);
RetCons.. Ret =g= DesRet;
Risk..	TotRisk =e= sum((I,J),X(I)*IndBeta(I,J));

Model MinVar /Fract,Pbeta,Pret,RetCons,Risk/;
Model MaxRet /Fract,Pbeta,Pret,BetaCons,Risk/;

if (ModelNum = 1 OR ModelNum = 3,
Solve MinVar using NLP minimizing TotRisk;
else
Solve MaxRet using NLP maximizing Ret;
);

scalar change /0/;
file progout    output file /portfolio.output/;
progout.pw=255;
progout.pc=6;
progout.nd=3;

put progout;

** finds whether current solution is feasible
change = 0;
Loop(I,
 if((X.l(I)<0),
  change = 3;
 );
);
if (sum(I,X.l(I))>1.0001 OR sum(I,X.l(I))<.9999,
 change=3;
);
if (ModelNum=2 or ModelNum=4,
Loop(J,
 if(DesBeta(J)+.001<Beta.l(J),
  change=3;
 );
);
elseif (ModelNum=1 or ModelNum=3),
 if (Ret.l<DesRet,
  change=3;
 );
);
** if change is 3, solution is infeasible

put / "*** Note that a lower or upper bound of 100 indicates an unbounded range ***";
put / ;
put / "                ","Factr","Bound","Ret","Risk";
put "Security proportions";
put / ;
if (change=0,
 put / "Optimal Solution";
elseif (change=3),
 put / "Problem Infeasible";
);
put " "," ";
put Ret.l,TotRisk.l;
Loop(I,
 put X.l(I);
);
put / ;

*********************************************************
* Sensitivity Analysis
* We want to do this on: DesBeta, DesRet, Factor


Alias (I,sI);
Alias (J,sJ);
Alias (J,lJ);
variable OptX(I) Temporary storage of optimal solution;
positive variable OptX;

Loop(sI,
 OptX.l(I) = X.l(I);
);


*********************************************************
* Changes in DesRet

scalar ULD /0/;
* specifies whether we are finding upper lim, lower lim, or done (0,1,2)
scalar delta	incremental change in diff /.1/;
scalar origret	original desired return;

origret = DesRet;

if (ModelNum = 1,
while (ULD < 2,
 change=0;
 delta = .1;
 DesRet = origret;

* Now increment by delta until either it appears to be unbounded 
* or the basis changes, or it becomes infeasible
 while (change=0,
  if (ULD = 0,
   DesRet = DesRet+delta;
  else
   DesRet = DesRet-delta;
  );
 Solve MinVar using NLP minimizing TotRisk;
* if DesRet is incrememted past the following ranges, assume unbounded
 if ((DesRet>1 or DesRet<0),
  change = 2;
 );
* check to see if new variable entered basis
 Loop(sI,
  if((OptX.l(sI)=0 AND X.l(sI)>0),
   change = 1;
  );
 );
* check for feasibility
 Loop(sI,
  if((X.l(sI)<0),
   change = 3;
  );
 );
 if (sum(sI,X.l(sI))>1.0001 OR sum(sI,X.l(sI))<.9999,
  change=3;
 );
 if(Ret.l<DesRet,
  change=3;
 );
 );

* If above loop find infeasibility or changed basis, 
* slowly change backwards to get close to change point
 if ((change=1 OR change=3),
 delta = -1*delta
* get delta within small range
  while ((delta<(-.01) OR delta>.01),
   if ((change = 1),
    delta = -1*abs(.5 * delta);
   elseif (change=3),
* change back to previous (feasible) value
    if (ULD = 0,
     DesRet = DesRet-delta;
    else
     DesRet = DesRet+delta;
    );
    delta = abs(.5 * delta);
   else
    delta = -.5 * delta;
   );
   if (ULD = 0,
    DesRet = DesRet+delta;
   else
    DesRet = DesRet-delta;
   );
   Solve MinVar using NLP minimizing TotRisk;
* now check for basis change or infeasibility
   change=0;
   Loop(sI,
    if((OptX.l(sI)=0 AND X.l(sI)>0),
     change = 1;
    );
   );

   Loop(sI,
    if((X.l(sI)<0),
     change = 3;
    );
   );
   if (sum(sI,X.l(sI))>1.0001 OR sum(sI,X.l(sI))<.9999,
    change=3;
   );
   if(Ret.l<DesRet,
    change=3;
   );

  );
* if it appears unbounded, output a value of 100
 elseif (change=2),
  DesRet = 100;
 );
* Output data to portfolio.output
   if (ULD=0,
    put / "Increase in DesRet";
   else
    put / "Decrease in DesRet";
   );
   put " ";
   put DesRet,Ret.l;
   put TotRisk.l;
   Loop(sI,
    put X.l(sI);
   );
 ULD = ULD + 1;
);
put / ;
);
* return DesRet to original value
DesRet = origret;

***********************************************************
* Changes in DesBeta

variable origbeta(sJ);

if (ModelNum = 2,
Loop(sJ,
origbeta.l(sJ) = DesBeta(sJ);
ULD = 0;
change = 0;

* find upper and lower bounds on DesBeta
while (ULD < 2,
 change = 0;
 DesBeta(sJ) = origbeta.l(sJ);
 delta =abs(.2 * DesBeta(sJ));

* Increment now until basis changes, becomes infeasible,
* or appears unbounded
 while (change=0,
  if (ULD = 0,
   DesBeta(sJ) = DesBeta(sJ)+delta;
  else
   DesBeta(sJ) = DesBeta(sJ)-delta;
  );
 Solve MaxRet using NLP maximizing Ret;
* check for unboundedness
 if ((DesBeta(sJ)>2 or DesBeta(sJ)<(-2)),
  change = 2;
 );
* check for new variable in basis
 Loop(sI,
  if((OptX.l(sI)=0 AND X.l(sI)>0),
   change = 1;
  );
 );
* check for infeasibility
 Loop(sI,
  if((X.l(sI)<0),
   change = 3;
  );
 );
 if (sum(sI,X.l(sI))>1.0001 OR sum(sI,X.l(sI))<.9999,
  change=3;
 );
 Loop(lJ,
  if(DesBeta(lJ)+.001<Beta.l(lJ),
   change=3;
  );
 );
 );

* if above loop found infeasibility or a new basis
* then slowly get within small range of change point
 if ((change=1 OR change=3),
 delta = -1*delta;
* while delta is outside of this range keep refining solution
  while ((delta<(-.02) OR delta>.02),
   if ((change = 1),
    delta = -1*abs(.5 * delta);
   elseif (change=3),
* change back to previous (feasible) value
    if (ULD = 0,
     DesBeta(sJ) = DesBeta(sJ)-delta;
    else
     DesBeta(sJ) = DesBeta(sJ)+delta;
    );
    delta = abs(.5 * delta);   
   else
    delta = -.5 * delta;
   );
   if (ULD = 0,
    DesBeta(sJ) = DesBeta(sJ)+delta;
   else
    DesBeta(sJ) = DesBeta(sJ)-delta;
   );
   Solve MaxRet using NLP maximizing Ret;
* now check for basis change or infeasibility
   change = 0;
   Loop(sI,
    If((OptX.l(sI)=0 AND X.l(sI)>0),
     change = 1;
    );
   );
  
   Loop(sI,
    if((X.l(sI)<0),
     change = 3;
    );
   );
   if (sum(sI,X.l(sI))>1.0001 OR sum(sI,X.l(sI))<.9999,
    change=3;
   );
   Loop(lJ,
    if(DesBeta(lJ)+.001<Beta.l(lJ),
     change=3;
    );
   );

  );
* if it appears unbounded, give DesBeta the value of 100
 elseif (change=2),
  DesBeta(sJ) = 100;
 );
* output data to portfolio.output
   if (ULD=0,
    put / "Increase in DesBeta ",sJ.tl;
   else
    put / "Decrease in DesBeta ",sJ.tl;
   );
   put DesBeta(sJ);
   put Ret.l,TotRisk.l;
   Loop(sI,
    put X.l(sI);
   );
 ULD = ULD + 1;
);
* return DesBeta to original value
DesBeta(sJ) = origbeta.l(sJ);
);
put / ;
);

****************************************************************
* Changes in factors


variable origfact(sJ);

if (ModelNum = 1 OR ModelNum = 2,
Loop(sJ,
origfact.l(sJ) = Factor(sJ);
ULD = 0;
change = 0;

while (ULD < 2,
 change = 0;
 delta = .02;
 Factor(sJ) = origfact.l(sJ);

* increment by delta until basis changes, becomes infeasible
* or problem appears unbounded
 while (change=0,
  if (ULD = 0,
   Factor(sJ) = Factor(sJ)+delta;
  else
   Factor(sJ) = Factor(sJ)-delta;
  );
  if (ModelNum = 1,
   Solve MinVar using NLP minimizing TotRisk;
  else
   Solve MaxRet using NLP maximizing Ret;
  );
* If outside of this range, assume Factor can be unbounded
 if ((Factor(sJ)>.2 or Factor(sJ)<(-.2)),
  change = 2;
 );
* Check for new basic variable
 Loop(sI,
  if((OptX.l(sI)=0 AND X.l(sI)>0),
   change = 1;
  );
 );
* check for infeasibility
 Loop(sI,
   if((X.l(sI)<0),
    change = 3;
   );
  );
 if (sum(sI,X.l(sI))>1.0001 OR sum(sI,X.l(sI))<.9999,
   change=3;
  );
 if (ModelNum=2,  
  Loop(lJ,
   if(DesBeta(lJ)+.001<Beta.l(lJ),
    change=3;
   );
  );
 elseif (ModelNum=1),
  if (Ret.l<DesRet,
   change=3;
  );
 );

 );

* if the above loop finds infeasibility or a new basis
* slowly refine solution to find change point
 if ((change=1 OR change=3),
 delta = -1*delta;
* stop when we are within to following range
  while ((delta<(-.005) OR delta>.005),
   if ((change = 1),
    delta = -1*abs(.5 * delta);
   elseif (change=3),
* change back to previous (feasible) value
    if (ULD = 0,
     Factor(sJ) = Factor(sJ)-delta;
    else
     Factor(sJ) = Factor(sJ)+delta;
    );
    delta = abs(.5 * delta);
   else
    delta = -.5 * delta;
   );
   if (ULD = 0,
    Factor(sJ) = Factor(sJ)+delta;
   else
    Factor(sJ) = Factor(sJ)-delta;
   );
   if (ModelNum = 1,
    Solve MinVar using NLP minimizing TotRisk;
   else
    Solve MaxRet using NLP maximizing Ret;
   );
* check for basis change or infeasibility
   change = 0;
   Loop(sI,
    If((OptX.l(sI)=0 AND X.l(sI)>0),
     change = 1;
    );
   );

   Loop(sI,
     if((X.l(sI)<0),
      change = 3;
     );
    );
    if (sum(sI,X.l(sI))>1.0001 OR sum(sI,X.l(sI))<.9999,
     change=3;
    );
   if (ModelNum=2,
    Loop(lJ,
     if(DesBeta(lJ)+.001<Beta.l(lJ),
      change=3;
     );
    );
   elseif (ModelNum=1),
    if (Ret.l<DesRet,
     change=3;
    );
   );

  );
* if problem appears unbounded, let Factor take value of 100
 elseif (change=2),
  Factor(sJ) = 100;
 );
* output data to portfolio.output
   if (ULD=0,
    put / "Increase in Factor ",sJ.tl;
   else
    put / "Decrease in Factor ",sJ.tl;
   );
   put Factor(sJ);
   put Ret.l,TotRisk.l;
   Loop(sI,
    put X.l(sI);
   );
 ULD = ULD + 1;
);
* return factor to original value
Factor(sJ) = origfact.l(sJ);
);
);


****************************************************************
putclose progout;