(!******************************************************
   Xpress-SP Example
   ======================
   ``````````````````
   FILE: StochasticStockBondAnalysis.mos
   
   DESCRIPTION: Stochastic multi-stage linear model
   
   TYPE: Finance
   
      We currently want to invest $b in stocks and bonds. 
   At the end of T-1 years, we would like to exceed the goal $G. 
   The investment is changed every year. The excess income would 
   imply an income of q% of the excess income and shortage would 
   imply borrowing at r% of the amount short.
   
   MODEL: 
		We also assume that at each stage the return on stocks and 
   bonds is 1.25 and 1.14 or 1.06 and 1.12 each with probability 0.5
   • Decision variables 
   	o x(i,t): amount invested in asset i at time t 
   	o y: shortage 
   	o w: excess 
   • Constraints 
   	o Balance of initial investment, returns and goal meeting requirements
   
   We also compare the results if all investments are done in stocks only.
   
   FURTHER INFO: see
   J.R.Birge and Francois Louveaux (97) - “Financial planning and control”
   - Introduction to Stochastic Programming, 20-24
      	
   DATE: Oct 2006
   
   (c) 2006 Dash Associates
       Author: Nitin Verma
       Dash Optimization Inc, NJ  
       
*******************************************************!)
model stochasticStockAndBondModel
uses 'mmsp'
parameters
	T=4
	Explicit=true
end-parameters
forward procedure generateTree
declarations
	Assets={"Stocks","Bonds"}
	Stages=1..T
	x:array(Assets,1..T-1) of spvar !amt invested
	w,y:spvar	!shortage, excess
	return:array(Assets,2..T) of sprand !return on investment
	b,G,q,r:real !initial wealth,goal(in$1000);excess,shortage(in%)
	ExpectedUtility:splinctr
	Balance:array(1..T) of splinctr
end-declarations

!data
b:=55;G:=80;q:=1;r:=4

spsetstages(Stages)

forall(i in Assets,t in 2..T) spsetstage(return(i,t),t)
generateTree

forall(i in Assets,t in 1..T-1) spsetstage(x(i,t),t)
spsetstage(y,T)
spsetstage(w,T)

ExpectedUtility:=q*y-r*w
forall(t in 1..T)
	if t=1 then Balance(t):=sum(i in Assets) x(i,t)=b;
	elif t=T then Balance(t):=sum(i in Assets) return(i,t)*x(i,t-1)-
							y+w=G;
	else Balance(t):=-sum(i in Assets) x(i,t)+
		sum(i in Assets) return(i,t)*x(i,t-1)=0;
	end-if
	
setparam("xsp_scen_based",false)
setparam('xsp_verbose',false)
maximize(ExpectedUtility)

ObjvalRP:=getobjval
writeln("Obj value in recourse problem = ",strfmt(ObjvalRP,0,2))
forall(t in 1..T-1) spfix(x("Bonds",t),0)
maximize(ExpectedUtility)
writeln("Obj value (all investments in Stocks only) = ",strfmt(getobjval,0,2))
writeln("VSS (if all investments in Stocks only) = ",strfmt(ObjvalRP-
		getobjval,0,2))
forall(t in 1..T-1) spunfix(x("Bonds",t))


procedure generateTree	
	declarations
		SetSrnd:set of sprand
		val:array(1..2,1..2) of real
		prob:array(1..2)  of real
		StockVal,BondVal,BranchProb:array(1..2) of real
	end-declarations
	
	if(Explicit) then 
		spcreatetree(2) !binary tree
		StockVal:=[1.25,1.06]
		BondVal:=[1.14,1.12]
		BranchProb:=[.5,.5]
		forall(t in 2..T) do !set Srnds
			spsetrand(return("Stocks",t),StockVal)
			spsetrand(return("Bonds",t),BondVal)
		end-do
		spsetprobcond(BranchProb) !set branch probabilities
		spgentree !generate tree
	else !use joint distribution functionality
		val:=[1.25, 1.06,
		  1.14, 1.12]
		prob:=[.5, .5]	
		forall(t in 2..T) do
			SetSrnd:={}
			forall(i in Assets) SetSrnd+={return(i,t)}	
			spsetdist(SetSrnd,val,prob) 
		end-do
		spgenexhtree
	end-if
		
end-procedure

end-model