MatLab code for qMu futures example

% Script prepared by Gethin Norman (sep03)
%
% script which returns maximum reward for investor
% requires file MINMAX.mat which is generated by PRISM and contains:
% S1, S2, S3 and S4 - transition matrices of the system 
% four matrices since at most four (non-deterministic) choices in a state
% (two for investor and two for environment)
% states - matrix of size n*6 representing the states of the system (n-number of states)
% (each row is of the form (i, v, p, c, b, m))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load files

load MINMAX;
fprintf('\nloaded matrices and states\n');

n=size(S1,1);
fprintf('\nnumber of states of the system: %d\n', n); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% construct matrices of investor and environment
% if m=0 then non-deterministic choices are made
% if m=1 then probabilistic choices are made

% NON_DETERMINISTIC CHOICES (when month variable is 0)
% start off with nothing (Cij i -choice of investor j choice of environment)
C11=sparse(n,n);
C12=sparse(n,n);
C21=sparse(n,n);
C22=sparse(n,n);

% when everything can choose - when the environment can bar the investor (b=0)
% the matrices have the following interpretation:
% S1 corresponds to investors first choice and environments first choice
% S2 corresponds to investors second choice and environments first choice
% S3 corresponds to investors first choice and environments second choice
% S4 corresponds to investors second choice and environments second choice
% fill rows with choices
C11(find(states(:,6)==0 & states(:,5)==0),:)=S1(find(states(:,6)==0 & states(:,5)==0),:);
C12(find(states(:,6)==0 & states(:,5)==0),:)=S3(find(states(:,6)==0 & states(:,5)==0),:);
C21(find(states(:,6)==0 & states(:,5)==0),:)=S2(find(states(:,6)==0 & states(:,5)==0),:);
C22(find(states(:,6)==0 & states(:,5)==0),:)=S4(find(states(:,6)==0 & states(:,5)==0),:);

% only investor can choose - when the environment cannot bar the investor (b=1)
% the matrices have the following interpretation:
% S1 corresponds to investors first choice 
% S2 corresponds to investors second choice 
% S3 and S4 empty
C11(find(states(:,6)==0 & states(:,5)==1),:)=S1(find(states(:,6)==0 & states(:,5)==1),:);
C12(find(states(:,6)==0 & states(:,5)==1),:)=S1(find(states(:,6)==0 & states(:,5)==1),:);
C21(find(states(:,6)==0 & states(:,5)==1),:)=S2(find(states(:,6)==0 & states(:,5)==1),:);
C22(find(states(:,6)==0 & states(:,5)==1),:)=S2(find(states(:,6)==0 & states(:,5)==1),:);

% PROBABILISTIC CHOICES (when month variable is 1)
% start off with nothing
P=sparse(n,n);

% the matrices have the following interpretation:
% S1 corresponds to probabilistic choices
% S2, S3 and S4 empty
% fill rows with choices
P(find(states(:,6)==1),:)=S1(find(states(:,6)==1),:);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fprintf('\nconstructed matrices of investor and environment\n');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CALCULATE COST
% cost in state where investor cashes in his/her shares (i=1 and m=0)
% at the same time remove outgoing transitions from these states

cost = sparse(n,1);

cost(find(states(:,1)==1 & states(:,6)==0 & states(:,5)==0),:)=states(find(states(:,1)==1 & states(:,6)==0 & states(:,5)==0),2);
C11(find(states(:,1)==1 & states(:,6)==0 & states(:,5)==0),:)=0;
C21(find(states(:,1)==1 & states(:,6)==0 & states(:,5)==0),:)=0;
C12(find(states(:,1)==1 & states(:,6)==0 & states(:,5)==0),:)=0;
C22(find(states(:,1)==1 & states(:,6)==0 & states(:,5)==0),:)=0;

fprintf('\ncalculated cost vectors and made target states absorbing\n');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ITERATIVE METHOD (min max)

% initial guess
res = sparse(n,1);
done = 0;
iter = 0;

fprintf('\nstarting iterations...');

% start iterations
while (done == 0)
	tmp = res;
	iter = iter + 1;
	tmp11 = cost+C11*P*res;
	tmp21 = cost+C21*P*res;
	tmp12 = cost+C12*P*res;
	tmp22 = cost+C22*P*res;
	% max for investor
	tmp1 = max(tmp11,tmp21);
	tmp2 = max(tmp12,tmp22);
	% min for environment
	res = min(tmp1,tmp2);
	% check convergence
	max(abs(res - tmp));
	if max(abs(res - tmp)) < 1*10^-8;
			done = 1;
	end
end

fprintf('finished after %d iterations (accuracy %d)\n', iter, 1*10^-8); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% these are the states we are interested (output from prism)
%262:(0,0,5,10,1,0) 
%746:(0,1,5,10,1,0) 
%1206:(0,2,5,10,1,0) 
%1622:(0,3,5,10,1,0) 
%1994:(0,4,5,10,1,0) 
%2322:(0,5,5,10,1,0) 
%2606:(0,6,5,10,1,0) 
%2846:(0,7,5,10,1,0) 
%3042:(0,8,5,10,1,0) 
%3194:(0,9,5,10,1,0) 
%3302:(0,10,5,10,1,0) 
% so need  to add one to each to get the corresponding position
% (since states indexed from 1 not 0 in matlab)

% get out results we require
res=res([263,747,1207,1623,1995,2323,2607,2847,3043,3195,3303],1);

fprintf('\nSOLUTION:', iter); 
for i=1:11
	fprintf('\nmaximum reward in state (i=0 & v=%d & p=0.5 & c=10 & b=1 & m=0) equals %1.6f', i-1, res(i)); 
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% check when max min

% iterative method
% initial guess
%res = sparse(n,1);
%done = 0;
%iter = 0;
%
%% max min
%while (done == 0)
%	tmp = res;
%	iter = iter + 1;
%
%	tmp11 = cost+C11*P*res;
%	tmp12 = cost+C21*P*res;
%	tmp12 = cost+C12*P*res;
%	tmp22 = cost+C22*P*res;
%
%	% max for environment
%	tmp1 = min(tmp11,tmp12);
%	tmp2 = min(tmp21,tmp22);
%	% min for environment
%	res = max(tmp1,tmp2);
%
%	% check convergence
%	max(abs(res - tmp))
%	if max(abs(res - tmp)) < 1*10^-8
%			done = 1;
%	end
%end