##
##                           dynamic.m
##
## Authors: Effrosyni Diamantoudi, Eiichi Miyagawa, and Licun Xue
## Date: Nov 15, 2007

## This program verifies that the strategy profile together with the belief
## system specified for one-sided commitment and patient workers is indeed
## sequentially rational. Verying the sequantional rationality manually is
## possible but extremely tedious. There are just many things to check. The
## following computer code checks the agents' incentives at every decision
## node and confirms that each agent has an incentive to follow the
## strategy. Thus by confirming the correctness of the code, you can verify
## the equilibrium indirectly. The code also generates a number of tables at
## the end of the process to summarize each agent's incentives. So you can
## also check the tables to see if the code is accurate.

## The computer code runs under Matlab and Octave. We have tested the code
## for Matlab 7.0.0 and Octave 2.9.4.

## Octave is a free reliable clone of Matlab and is available for many
## computer systems from http://www.octave.org/. For users of Microsoft
## Windows, a nice interface is provided by Octave Workshop, which is also
## freely available at http://www.math.mcgill.ca/loisel/octave-workshop/.

## To run the code, start Octave or Matlab. Put this file in the current folder
## and type "dynamic". The code finishes within a few seconds and gives messages
## confirming the incentives. The code also generates two files named
## "TablesForFirms.txt" and "TablesForWorkers.txt", which summarize the incentives
## of firms and workers, respectively.

## Please send us any question or comment by emailing to miyagawa@econ.kobe-u.ac.jp
##############################################################################
clear;


### Static Physical Data

F = 3; # Number of Firms
W = 3; # Number of Workers
f = [1 2 3]'; # f(1) = agent 1, f(2) = agent 2, f(3) = agent 3
w = [4 5 6]'; # w(1) = agent 4, w(2) = agent 5, w(3) = agent 6



### Preference Data

delta = 0.9; # discount factor
# The strategy profile is an equilibrium if delta > (0.7)^(1/3) = 0.887905.

# Period-payoffs
# PeriodPayoff(i,j) = i's period-utility from agent j.
# We normalize PeriodPayoff(i,i) = 0.
PeriodPayoff = -100*ones(F+W) + 100*eye(F+W);
# u(k) = utilify from the kth choice (k = 1, 2, 3).
# This applies to every agent.
u = [100 70 40];
#                        w1   w2   w3
PeriodPayoff(f(1),w') = [u(1) u(3) u(2)]; # firm 1's period-payoff function
PeriodPayoff(f(2),w') = [u(3) u(2) u(1)]; # firm 2
PeriodPayoff(f(3),w') = [u(2) u(3) u(1)]; # firm 3
#                        f1   f2   f3
PeriodPayoff(w(1),f') = [u(3) u(1) u(2)]; # worker 1
PeriodPayoff(w(2),f') = [u(3) u(2) u(1)]; # worker 2
PeriodPayoff(w(3),f') = [u(1) u(2) u(3)]; # worker 3



### States

S = 34;		# Number of States
Matching = [1:F+W]'*ones(1,S);
# Matching(i,s) = i's partner in state s
Matching(w(1),:) = [#
    w(1); f(1); w(1); w(1); f(2);
	f(3); f(3); f(3); w(1); f(1);
	f(1); f(1); w(1); w(1); w(1);
	f(2); f(2); f(2); w(1); w(1);
	f(1); w(1); f(3); f(3); f(3);
	w(1); w(1); w(1); f(3); f(2);
	f(2); f(2); f(1); f(1)      ]';
Matching(w(2),:) = [#
	w(2); w(2); f(3); f(3); f(3);
	w(2); w(2); f(2); f(2); w(2);
	f(2); f(3); w(2); w(2); w(2);
	w(2); f(3); w(2); f(3); f(2);
	w(2); f(2); f(2); f(1); f(1);
	f(1); f(1); f(1); w(2); f(1);
	w(2); f(1); f(3); f(2)     ]';
Matching(w(3),:) = [#
	w(3); f(2); f(2); f(1); f(1);
	f(2); f(1); f(1); f(3); f(3);
	w(3); w(3); f(3); f(2); f(1);
	f(1); w(3); w(3); w(3); f(1);
	w(3); w(3); w(3); w(3); f(2);
	w(3); f(2); f(3); w(3); w(3);
	f(3); f(3); f(2); f(3)      ]';

# Find each firm's partner by using mu(mu(i)) = i
for state = 1:S
	for j = 1:W
		Matching(Matching(w(j),state),state) = w(j);
	end;
end;

# PeriodPayoffInState(s,i) = Agent i's period-payoff in state s
PeriodPayoffInState = zeros(S,F+W);
for state = 1:S
	for i = 1:F+W
	PeriodPayoffInState(state,i) = PeriodPayoff(i,Matching(i,state));
	end;
end;



### Equilibrium

## Each firm's offer in each state
## EqmOffer(f(i), s) = f(i)'s equilibrium offer in state s.
EqmOffer = zeros(F,S);
EqmOffer(f(1),:) = [#
	w(1); w(1); w(3); w(3); w(3);
	w(3); w(3); w(3); w(1); w(1);
	w(1); w(1); w(1); w(3); w(3);
	w(3); w(3); w(3); w(3); w(3);
	w(1); w(1); w(3); w(2); w(2);
	w(2); w(2); w(2); w(3); w(2);
	w(3); w(2); w(1); w(1)			]';
EqmOffer(f(2),:) = [#
	w(2); w(3); w(3); w(1); w(1);
	w(3); w(2); w(2); w(2); w(2);
	w(2); w(3); w(2); w(3); w(1);
	w(1); w(1); w(1); w(1); w(2);
	w(2); w(2); w(2); w(2); w(3);
	w(2); w(3); w(2); w(2); w(1);
	w(1); w(1); w(3); w(2)			]';
EqmOffer(f(3),:) = [#
	w(3); w(2); w(2); w(2); w(2);
	w(1); w(1); w(1); w(3); w(3);
	w(3); w(2); w(3); w(2); w(2);
	w(2); w(2); w(2); w(2); w(1);
	w(3); w(3); w(1); w(1); w(1);
	w(3); w(2); w(3); w(1); w(3);
	w(3); w(3); w(2); w(3)			]';


## Each worker's cutoff
# If Cutoff(w(i), state) = k, it means that w(i)'s action in the state
# is to accept the best offer that is at least as good as agent k.
Cutoff = [1:F+W]'*ones(1,S);

# The default cutoff level is oneself.
# Firms are included in the rows just for convenience. This part will not be used.

# When the cutoff level > oneself
Cutoff(w(1), [6 12 14 15 19 27]) = f(2);
Cutoff(w(1), [2 3 33]) = f(3);
Cutoff(w(2), [31 32]) = f(2);
Cutoff(w(2), [15 18 30]) = f(3);
Cutoff(w(3), [24 25]) = f(1);
Cutoff(w(3), [14]) = f(2);

# When the # of offers > 1, w(1) changes the cutoff in states 6 and 14
CutoffSpecial = Cutoff;
CutoffSpecial(w(1),[6 14]) = f(3);

# The cutoff in terms of period-utilities
CutoffUtilityLevel = zeros(2*W,S);
for state=1:S
for k=1:W
	CutoffUtilityLevel(k,state)   = PeriodPayoff(w(k),Cutoff(w(k),state));
	CutoffUtilityLevel(W+k,state) = PeriodPayoff(w(k),CutoffSpecial(w(k),state));
end;
end;


## Workers' belief when the expected offer did not arrive.
#  At every such situation, the worker believes that
#  the firm did not make any offer to any worker.
#  This belief is used when PotentialOffer(i) is computed below.


## State Transition

# NextState(s) = the state that comes after state s in equilibrium.
NextState = [#
34 3 4 5 5 15 8 8 34 34 34 3 34 4 5 5 5 5 5 8 34 34 8 23 24 9 3 9 8 31 16 31 3 34];

# Transition matrix
# StateTransition(s,t) = 1 if state t comes after state s in equilibrium,
# StateTransition(s,t) = 0 otherwise.
StateTransition = zeros(S,S);
for state = 1:S
	StateTransition(state, NextState(state)) = 1;
end;

# Check if the inverse exists
if (det(eye(S) - delta*StateTransition) == 0); disp('Determinant is zero'); end;

# Continuation value from each state
# Value(state, i) = agent i's continuation value from the state.
Value = (1-delta)*inv(eye(S) - delta*StateTransition)*PeriodPayoffInState;



###########################################
### Check the optimality of firms' strategy
###########################################

## Initialize before iterations

IsEqmActionIndeedOptimal = -1.*ones(F+W,S);  # The initial value = -1
ContingentStateForFirms = zeros(F*(W+1),S);

## Check each state
for state = 1:S;


## The matching at the beginning of the period
mu = Matching(:,state);

InactiveFirms = find(mu(f) ~= f)';
ActiveFirms   = find(mu(f) == f)';

## Renewal offers from inactive firms
RenewalOffers = -1.*ones(F, F+W);
for i = InactiveFirms
	RenewalOffers(f(i), mu(f(i))) = 1;
end;

## Choose any firm
for i = ActiveFirms

## Registor the other firms' equilibrium offer
OthersOffers = RenewalOffers;
for k = ActiveFirms
	if (k ~= i)
		OthersOffers(f(k), EqmOffer(f(k),state)) = 1;
	end;
end;

ContingentUtility = zeros(F+W, 1);

## Choose any possible offer
PossibleOffers = [w' i]; ## The set of offers that the firm can make
for offer = PossibleOffers

## All the outstanding offers from the firms
AllOffers = OthersOffers;
AllOffers(f(i), offer) = 1;

## Options for each worker (in terms of period-utilities)
UtilityOptionsForWorkers = -1.*ones(F+W,W);
for h = 1:W
	UtilityOptionsForWorkers(f,h) = AllOffers(f, w(h)).*PeriodPayoff(w(h),f)';
	UtilityOptionsForWorkers(w(h),h) = 0; # The option of being unemployed
end;

## Best offer for each worker
[Maxima, Maximizer] = max(UtilityOptionsForWorkers);
## Maxima(i) = the maximum period-utility for worker i given the offers.
## Maximizer(i) = the most preferred offer for worker i.

## New matching after the workers' choice
NewMatching = [1: F+W]';
for h = 1:W
	special = (h==1)*( (state==6) | (state==14) )*( sum( AllOffers(:,w(h)) == 1 ) > 1 );
		# if it is worker 1 in state 6 or 14,
		# and if the worker has more than one offer,
		# then special = 1, which makes the worker's cutoff level equal CutoffSpecial.
	if (Maxima(h) < CutoffUtilityLevel(special*W + h, state))
		NewMatching(w(h)) = w(h);  ## The best offer does not meet the cutoff
	elseif (Maxima(h) == 0)		   ## Maxima could be 0 if there was no offer
		NewMatching(w(h)) = w(h);
	else
		NewMatching(w(h)) = Maximizer(h);  ## The cutoff is met.
	end;
end;
# Compute each firm's partner by using mu(mu(i)) = i.
for h = 1:W
	NewMatching(NewMatching(w(h))) = w(h);
end;

## The new state number (1,...,34)
NewStateCandidates = any( Matching - NewMatching*ones(1,S) ) == 0;
	# NewStateCandidates(s) = 1 if NewMatching coincides with Matching(:,s),
	#                         0 otherwise.
if (sum(NewStateCandidates) ~= 1); disp('Problem with NewStateCandidates.'); end;
NewState = NewStateCandidates*[1:S]'; # The unique s such that NewStateCandidates(s)=1

## Record the contingent state

IndexOfOffer = find( PossibleOffers == offer );
ContingentStateForFirms( (i-1)*(W+1) + IndexOfOffer, state ) = NewState;

## By computing the value from the new state,
## we can compute the value from making the particular offer.
ContingentUtility(offer) = Value(NewState,f(i));

end; # offer

[HighestUtility, BestOfferToMake] = max(ContingentUtility);

if (BestOfferToMake ~= EqmOffer(f(i), state))
	IsEqmActionIndeedOptimal(f(i),state) = 0;
	error('suboptimal');
	# The best offer to make differs from the prescribed action.
else
	IsEqmActionIndeedOptimal(f(i),state) = 1;
	# This is what should be.
end;

end; # i
end; # state

## Report whether the strategy is indeed optimal.
if sum( IsEqmActionIndeedOptimal == 0 ) > 0
	disp('Oops... firms are not behaving optimally!');
else
	disp('Firms are indeed behaving optimally.');
end;

StatesWhereFirmsCanMove = zeros(F,S);
StatesWhereFirmsCanMove(1,:) = any( ContingentStateForFirms([1:W+1],:) );
StatesWhereFirmsCanMove(2,:) = any( ContingentStateForFirms( W+1 + [1:W+1], :) );
StatesWhereFirmsCanMove(3,:) = any( ContingentStateForFirms( 2*(W+1) + [1:W+1], :) );



#############################################
### Check the optimality of workers' strategy
#############################################

PossibleOfferProfiles = [#
	0 0 0;
	1 0 0;
	0 1 0;
	0 0 1;
	1 1 0;
	1 0 1;
	0 1 1;
	1 1 1 ];
NumberOfOfferProfiles = size(PossibleOfferProfiles,1);
OptimalContingentChoice = zeros(W*S, NumberOfOfferProfiles);
ContingentStateForWorkers = zeros(W*S*(F+1), NumberOfOfferProfiles);


## Choose any state
for state = 1:S

mu = Matching(:,state); 	# The matching at the beginning of the state
PrescribedOffer = EqmOffer(:, state); 	# For each agent
OneIfActiveFirm = (mu(f) == f)';		# OneIfActiveFirm(i) = 1 if f(i) is active
InactiveFirms 	= find( mu(f) ~= f )';	# List of inactive firms


## Renewal offers from inactive firms
RenewalOffers = -1.*ones(F, F+W);
for i = InactiveFirms
	RenewalOffers(f(i), mu(f(i))) = 1;
end;

## Choose any worker
for j = 1:W

PotentialOffer = zeros(1,F);
AllOffers = RenewalOffers;

for OfferByFirm1 = 0:OneIfActiveFirm(1)
# Each active firm has 2 possibilities that matter for worker j:
# equilibrium action and the relevant deviation

PotentialOffer(1) = ( 1-OneIfActiveFirm(1) )*mu(f(1)) + OneIfActiveFirm(1)*(1-OfferByFirm1)*PrescribedOffer(1) + OneIfActiveFirm(1)*OfferByFirm1*[ ( PrescribedOffer(1)==w(j) )*f(1) + ( PrescribedOffer(1) ~= w(j) )*w(j) ];
# Firm 1 makes an offer to:
#	mu(f(1)) if the firm is inactive;
# 	PrescribedOffer(1) if the firm is active and does not deviate;
#	itself if the firm is prescribed to make an offer to w(j) but deviates
#	       and makes no offer to any worker;
#   (Worker j believes this way if she does not receive the offer from the firm.)
#	w(j)   if the firm is prescribed to make an offer elsewhere but deviates
#          and makes an offer to worker j.

AllOffers(f(1), :) = RenewalOffers(f(1),:);
AllOffers(f(1), PotentialOffer(1)) = 1;

# The same possibilities for f(2)
for OfferByFirm2 = 0:OneIfActiveFirm(2)

PotentialOffer(2) = ( 1-OneIfActiveFirm(2) )*mu(f(2)) + OneIfActiveFirm(2)*(1-OfferByFirm2)*PrescribedOffer(2) + OneIfActiveFirm(2)*OfferByFirm2*[ ( PrescribedOffer(2)==w(j) )*f(2) + ( PrescribedOffer(2) ~= w(j) )*w(j) ];

AllOffers(f(2), :) = RenewalOffers(f(2),:);
AllOffers(f(2), PotentialOffer(2)) = 1;

# The same possibilities for f(3)
for OfferByFirm3 = 0:OneIfActiveFirm(3)

PotentialOffer(3) = ( 1-OneIfActiveFirm(3) )*mu(f(3)) + OneIfActiveFirm(3)*(1-OfferByFirm3)*PrescribedOffer(3) + OneIfActiveFirm(3)*OfferByFirm3*[ ( PrescribedOffer(3)==w(j) )*f(3) + ( PrescribedOffer(3) ~= w(j) )*w(j) ];

AllOffers(f(3), :) = RenewalOffers(f(3),:);
AllOffers(f(3), PotentialOffer(3)) = 1;


## Options for each worker (in terms of period-utilities)
UtilityOptionsForWorkers = -1.*ones(F+W,W);
for h = 1:W
	UtilityOptionsForWorkers(f,h) = AllOffers(f, w(h)).*PeriodPayoff(w(h),f)';
	UtilityOptionsForWorkers(w(h),h) = 0; # The option of being unemployed
end;

## Best offer for each worker
[Maxima, Maximizer] = max(UtilityOptionsForWorkers);

## New matching after the workers' choice
NewMatching = [1: F+W]';
for h = 1:W
	special = (h==1)*( (state==6) | (state==14) )*( sum( AllOffers(:,w(h)) == 1 ) > 1 );
		# if it is worker 1 in state 6 or 14,
		# and if the worker has more than one offer,
		# then the worker's cutoff level is CutoffSpecial.
	if (   Maxima(h) < CutoffUtilityLevel(special*W + h, state)    )
		NewMatching(w(h)) = w(h);  ## The best offer does not meet the cutoff
	elseif (Maxima(h) == 0)		   ## Maxima could be 0 if there was no offer
		NewMatching(w(h)) = w(h);
	else
		NewMatching(w(h)) = Maximizer(h);  ## The cutoff is met.
	end;
end;

OfferToTakeInEqm = NewMatching(w(j));	# The offer that w(j) takes in equilibrium.
NewMatchingForOthers = NewMatching; 	# The other workers' eqm actions are recorded.

# The set of offers available to w(j)
OffersForWorkerJ = [ find( AllOffers(:, w(j))' == 1) w(j) ];
# Initialize
ContingentUtilityForWorkerJ = -1.*ones(1,F+W);

## Choose any offer that worker j may choose
for offer = OffersForWorkerJ

## The resulting matching
NewMatching = NewMatchingForOthers;
NewMatching(w(j)) = offer;
for h = 1:W
	NewMatching(NewMatching(w(h))) = w(h); ## if w(h) gets f(k), then f(k) gets w(h)
end;

## The new state number (1,...,34)
NewStateCandidates = any( Matching - NewMatching*ones(1,S) ) == 0;
	# NewStateCandidates(s) = 1 if NewMatching coincides with Matching(:,s),
	#                         0 otherwise.
if (sum(NewStateCandidates) ~= 1); disp('Problem with NewStateCandidates.'); end;
NewState = NewStateCandidates*[1:S]'; # The unique s such that NewStateCandidates(s)=1

## By computing the value from the new state,
## we compute the continuation value from taking the particular offer
ContingentUtilityForWorkerJ(offer) = Value(NewState,w(j));

## Record-keeping
OfferSetNumber = find( any( PossibleOfferProfiles' - (AllOffers(:,w(j)) == 1)*ones(1,NumberOfOfferProfiles) ) == 0);

ContingentStateForWorkers((j-1)*S*(F+1) + (state-1)*(F+1) + min(offer,F+1),OfferSetNumber)=NewState;

end; # offer

## Identify the best action for w(j).
[HighestUtility, BestOfferToTake] = max(ContingentUtilityForWorkerJ);

## Stop the program if the best action differs from the prescribed action.
if (BestOfferToTake ~= OfferToTakeInEqm)
	IsEqmActionIndeedOptimal(w(j),state) = 0;
	error('suboptimal');
else
	IsEqmActionIndeedOptimal(w(j),state) = abs(IsEqmActionIndeedOptimal(w(j),state));
end;

## Record the best action for later use.
OptimalContingentChoice( (j-1)*S + state, OfferSetNumber) = min( BestOfferToTake, F+1);

end; end; end;
end; # j
end; # state

## Report whether each worker's strategy is indeed optimal.
if sum( IsEqmActionIndeedOptimal(w,:) == 0 ) > 0
	disp('Oops... Workers are not behaving optimally!');
else
	disp('Workers are indeed behaving optimally.');
end;



#######################################################
### Print the firms' incentives to "TablesForFirms.txt"
#######################################################


## Open a file for recording
fp = fopen('TablesForFirms.txt','wt');

## Print a note.
fprintf(fp,'How to read the tables:\n1. xx ( yy) in a cell says if the firm makes an offer to the worker\n   of the column in the state of the row, state xx is induced,\n   from which the continuation value to the firm is yy.\n2. * denotes the firm''s equilibrium offer in the state.\n3. The tables omit the states where the firm cannot move.\n');

## Pick each firm
for i = 1:F

fprintf(fp,'\n');
fprintf(fp,'Firm %1.0f\n', i);
fprintf(fp,'State#  w(1)       w(2)       w(3)       none\n');

## Pick each state
for state = 1:S

## Check if firm i can move in the state
if ( StatesWhereFirmsCanMove(i,state) == 1 )

## The contingent state for each offer the firm may make.
ContState = ContingentStateForFirms((i-1)*(W+1) + [1:W+1], state);

## Put a * on the equilibrium offer
eqm = EqmOffer(i,state) - F;
if (eqm <= 0); eqm = 4; end;
star = [' ';' ';' ';' '];
star(eqm,:) = '*';

## Print out
fprintf(fp,'%-7.0f %2.0f (%3.0f)%s  %2.0f (%3.0f)%s  %2.0f (%3.0f)%s  %2.0f (%3.0f)%s\n', state, ContState(1),Value(ContState(1),i),star(1,:),ContState(2),Value(ContState(2),i),star(2,:),ContState(3),Value(ContState(3),i),star(3,:),ContState(4),Value(ContState(4),i),star(4,:));

end; # if (...)
end; # state
end; # i

## Close the file.
fclose(fp);



#######################################################
### Print the workers' incentives to "TablesForWorkers"
#######################################################

## Open a file for recording.
fw = fopen('TablesForWorkers.txt','wt');

## Print a note.
fprintf(fw,'How to read the tables:\n1. xx ( yy) in a cell means that if the worker receives the set of offers\n   denoted by the row and chooses the offer denoted by the column, then\n   it induces state xx, from which the continuation value to the worker is yy.\n2. * denotes the equilibrium choice in the contingency of the row.\n');

## Labels for rows
OfferSet = [#
	'none   ';
	'{f1}   ';
	'{f2}   ';
	'{f3}   ';
	'{f1,f2}';
	'{f1,f3}';
	'{f2,f3}';
	'all    '];

## Labels for columns
OfferName = [#
	'f1';
	'f2';
	'f3';
	'no'];

## Choose each state
for state=1:S

fprintf(fw,'\n');
fprintf(fw,'State %1.0f\n', state);
fprintf(fw,'           f1         f2         f3         none\n');

## Choose each worker
for j=1:W
fprintf(fw,'Worker %1.0f\n', j);

## Choose each possible offer set (Ignore the empty set)
for h = 2:8
fprintf(fw,'%-8s:  ',OfferSet(h,:));

## The contingent state for each choice of the worker: 4-dim vector
ContState = ContingentStateForWorkers((j-1)*S*(F+1) + (state-1)*(F+1) + [1:F+1], h);

## The optimal choice for the worker
OptimalChoice = OptimalContingentChoice( (j-1)*S + state,h);

## Put a star on the optimal choice
Star = [' ';' ';' ';' '];
if (OptimalChoice > 0); Star(OptimalChoice,:) = '*'; end;

## Choose each option (a firm or reject all)
for k = 1:4

## Print "--" if the option is not available for the worker.
if (ContState(k) == 0); fprintf(fw,'--         ');

## Otherwise, print the state that the option induces.
else
fprintf(fw,'%2.0f (%3.0f)%s  ',ContState(k),Value(ContState(k),w(j)),Star(k,:));
end;

end; # k
fprintf(fw,'\n');
end; # h
fprintf(fw,'\n');
end; # j
end; # state

## Close the file.
fclose(fw);


### Final Message
disp('');
disp('See TablesForFirms.txt and TablesForWorkers.txt, which have been');
disp('generated in the current folder and summarize the incentives');
disp('of firms and workers, respectively.');
disp('');