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