Numerics for a game theory calculation using expected utility

Question

I am trying to replicate Bruce B. de Mesquita's (BDM) results on political game theory for prediction. Based on where actors stand on issues, their capabilities, salience, BDM's method attempts to find the eventual decision point by simulating a game. He reportedly used this method with much success; and published his results in successive journals, the latest of which is (1). This is his so-called "expected utility method", there is a newer method (3) but there is less documentation on that, so I wanted to use EU model first.

Scholz et.al tried to replicate the findings and documented his work here (2). I took his work as basis, since a lot of BDM articles / books are behind paywalls. There are also the gentleman here (4), they took Scholz's work as the basis, added a machine learning method on top, and created a new product.

I wrote the code, however I am not sure I was successful at replicating results.

import pandas as pd
import numpy as np
import itertools

Q = 1.0 ; T = 1.0

class Game:

    def __init__(self,df):
        self.df = df
        self.df_capability = df.Capability.to_dict()    
        self.df_position = df.Position.to_dict()    
        self.df_salience = df.Salience.to_dict()    
        self.max_pos = df.Position.max()
        self.min_pos = df.Position.min()

    def weighted_median(self):
        self.df['w'] = self.df.Capability*self.df.Salience 
        self.df['w'] = self.df['w'] / self.df['w'].sum()
        self.df['w'] = self.df['w'].cumsum()
        return float(self.df[self.df['w']>=0.5].head(1).Position)

    def mean(self):
        return (self.df.Capability*self.df.Position*self.df.Salience).sum() / \
               (self.df.Capability*self.df.Salience).sum()

    def Usi_i(self,i,j,ri=1.):
        tmp1 = self.df_position[i]-self.df_position[j]
        tmp2 = self.max_pos-self.min_pos
        return 2. - 4.0 * ( (0.5-0.5*np.abs(float(tmp1)/tmp2) )**ri)

    def Ufi_i(self,i,j,ri=1.):
        tmp1 = self.df_position[i]-self.df_position[j]
        tmp2 = self.df.Position.max()-self.df.Position.min()
        return 2. - 4.0 * ( (0.5+0.5*np.abs(float(tmp1)/tmp2) )**ri )

    def Usq_i(self,i,ri=1.):
        return 2.-(4.*(0.5**ri))

    def Ui_ij(self,i,j):
        tmp1 = self.df_position[i] - self.df_position[j]
        tmp2 = self.max_pos-self.min_pos
        return 1. - 2.*np.abs(float(tmp1) / tmp2) 

    def v(self,i,j,k):
        return self.df_capability[i]*self.df_salience[i]*(self.Ui_ij(i,j)-self.Ui_ij(i,k)) 

    def Pi(self,i):
        l = np.array([[i,j,k] for (j,k) in itertools.combinations(range(len(self.df)), 2 ) if i!=j and i!=k])
        U_filter = np.array(map(lambda (i,j,k): self.Ui_ij(j,i)>self.Ui_ij(i,k), l))
        lpos = l[U_filter]
        tmp1 = np.sum(map(lambda (i,j,k): self.v(j,i,k), lpos))
        tmp2 = np.sum(map(lambda (i,j,k): self.v(j,i,k), l))
        return float(tmp1)/tmp2

    def Ubi_i(self,i,j,ri=1):
        tmp1 = np.abs(self.df_position[i] - self.weighted_median()) + \
               np.abs(self.df_position[i] - self.df_position[j])
        tmp2 = np.abs(self.max_pos-self.min_pos)
        return 2. - (4. * (0.5 - (0.25 * float(tmp1) / tmp2))**ri)

    def Uwi_i(self,i,j,ri=1):
        tmp1 = np.abs(self.df_position[i] - self.weighted_median()) + \
               np.abs(self.df_position[i] - self.df_position[j])
        tmp2 = np.abs(self.max_pos-self.min_pos)
        return 2. - (4. * (0.5 + (0.25 * float(tmp1) / tmp2))**ri)

    def EU_i(self,i,j,r=1):
        term1 = self.df_salience[j] * \
                ( self.Pi(i)*self.Usi_i(i,j,r) + ( 1.-self.Pi(i) )*self.Ufi_i(i,j,r) )
        term2 = (1-self.df_salience[j])*self.Usi_i(i,j,r)
        #term3 = -self.Qij(j,i)*self.Usq_i(i,r)
        #term4 = -(1.-self.Qij(j,i))*( T*self.Ubi_i(i,j,r) + (1.-T)*self.Uwi_i(i,j,r) )
        term3 = -Q*self.Usq_i(i,r)
        term4 = -(1.-Q)*( T*self.Ubi_i(i,j,r) + (1.-T)*self.Uwi_i(i,j,r) )
        return term1+term2+term3+term4

    d ef EU_j(self,i,j,r=1):
        return self.EU_i(j,i,r)

    def Ri(self,i):
        # get all j's except i
        l = [x for x in range(len(self.df)) if x!= i]
        tmp = np.array(map(lambda x: self.EU_j(i,x), l))
        numterm1 = 2*np.sum(tmp)
        numterm2 = (len(self.df)-1)*np.max(tmp)
        numterm3 = (len(self.df)-1)*np.min(tmp)
        return float(numterm1-numterm2-numterm3) / (numterm2-numterm3)

    def ri(self,i):
        Ri_tmp = self.Ri(i)
        return (1-Ri_tmp/3.) / (1+Ri_tmp/3.)

    def Qij(self,i,j):
        l = np.array([k for k in range(len(self.df))])
        res = map(lambda x: self.Pi(k)+(1-self.df_salience[k]),l)
        return np.product(res)

    def do_round(self,df):
        self.df = df; df_new = self.df.copy()        
        # reinit
        self.df_capability = self.df.Capability.to_dict()    
        self.df_position = self.df.Position.to_dict()    
        self.df_salience = self.df.Salience.to_dict()    
        self.max_pos = self.df.Position.max()
        self.min_pos = self.df.Position.min()

        offers = [list() for i in range(len(self.df))]
        ris = [self.ri(i) for i in range(len(self.df))]
        for (i,j) in itertools.combinations(range(len(self.df)), 2 ):
            eui = self.EU_i(i,j,r=ris[i])
            euj = self.EU_j(i,j,r=ris[j])
            if eui > 0 and euj > 0:
                # conflict
                mid_step = (self.df_position[i]-self.df_position[j])/2.
                print i,j,eui,euj,'conflict, both step', mid_step, -mid_step
                offers[j].append(mid_step)
                offers[i].append(-mid_step)
            elif eui > 0 and euj < 0 and np.abs(eui) > np.abs(euj):
                # compromise - actor i has the upper hand
                print i,j,eui,euj,'compromise', i, 'upper hand'
                xhat = (self.df_position[i]-self.df_position[j]) * np.abs(euj/eui)
                offers[j].append(xhat)
            elif eui < 0 and euj > 0 and np.abs(eui) < np.abs(euj):
                # compromise - actor j has the upper hand
                print i,j,eui,euj,'compromise', j, 'upper hand'
                xhat = (self.df_position[j]-self.df_position[i]) * np.abs(eui/euj)
                offers[i].append(xhat)
            elif eui > 0 and euj < 0 and np.abs(eui) < np.abs(euj):
                # capinulation - actor i has upper hand
                j_moves = self.df_position[i]-self.df_position[j]
                print i,j,eui,euj,'capitulate', i, 'wins', j, 'moves',j_moves
                offers[j].append(j_moves)
            elif eui < 0 and euj > 0 and np.abs(eui) > np.abs(euj):
                # capitulation - actor j has upper hand
                i_moves = self.df_position[j]-self.df_position[i]
                print i,j,eui,euj,'capitulate', j, 'wins', i, 'moves',i_moves
                offers[i].append(i_moves)
            else:
                print i,j,eui,euj,'nothing'

        print offers
        df_new['offer'] = map(lambda x: 0 if len(x)==0 else x[np.argmin(np.abs(x))],offers)
        df_new.loc[:,'Position'] = df_new.Position + df_new.offer
        df_new.loc[df_new['Position']>self.max_pos,'Position'] = self.max_pos
        df_new.loc[df_new['Position']<self.min_pos,'Position'] = self.min_pos
        return df_new

To run, there is run.py:

import pandas as pd, sys
import numpy as np, matplotlib.pylab as plt
import scholz, itertools

if len(sys.argv) < 3:
    print "\nUsage: run.py [CSV] [ROUNDS]"
    exit()

df = pd.read_csv(sys.argv[1]); print df
df.Position = df.Position.astype(float)
df.Capability = df.Capability.astype(float)
df.Salience = df.Salience/100.

game = scholz.Game(df)

results = pd.DataFrame(index=df.index)
for i in range(int(sys.argv[2])):
    results[i] = df.Position
    df = game.do_round(df)
    print df
    print 'weighted_median', game.weighted_median(), 'mean', game.mean()

results =  results.T
results.columns = df.Actor
print results
results.plot()
plt.savefig('out-%s.png' % sys.argv[1])

I ran this code on EU emission agreement, Iran presidential election data from (4), on the British EMU data from (5) (for Labor party case), and two small synthetic datasets I created.

Actor,Capability,Position,Salience
Netherlands,8,40,80
Belgium,8,70,40
Luxembourg,3,40,20
Germany,16,40,80
France,16,100,60
Italy,16,100,60
UK,16,100,90
Ireland,5,70,10
Denmark,5,40,100
Greece,8,70,70

Actor,Capability,Position,Salience
Jalili,24,10,70
Haddad,8,20,100
Gharazi,1,40,100
Rezayi,20,40,60
Ghalibaf,64,50,100
Velayati,7,50,25
Ruhani,21,80,100
Aref,30,100,70

Actor,Capability,Position,Salience
Labor Pro EMU,100,75,40
Labor Eurosceptic,50,35,40
The Bank of England,10,50,60
Technocrats,10,95,40
British Industry,10,50,40
Institute of Directors,10,40,40
Financial Investors,10,85,60
Conservative Eurosceptics,30,5,95
Conservative Europhiles,30,60,50

Actor,Capability,Position,Salience
A,100,100,100
B,100,90,100
C,50,50,50
D,5,5,10
E,10,10,20

Actor,Capability,Position,Salience
A,100,5,100
B,100,10,100
C,50,50,50
D,5,100,10
E,10,90,20

For EU emission (4) reports the result should have been around 8, I get 6.5. For Iran outcome is around 60, favoring reformers but this is far cry from Preana's and BDMs findings which is around 80. For EMU data, authors report anti-euro finding near 4, my finding is around 60.

The synthetic dataset is fine, always coalescing near top and bottom, but this is a simple case. I am attaching the graph outputs below as well.

---

1. Bueno De Mesquita BB (1994) Political forecasting: an expected utility method. In: Stockman F (ed.) European Community Decision Making. Yale, CT: Yale University Press, Chapter 4, 71–104.
2. https://oficiodesociologo.files.wordpress.com/2012/03/scholz-et-all-unravelling-bueno-de-mesquita-s-group-decision-model.pdf
3. A New Model for Predicting Policy Choices: Preliminary Tests http://irworkshop.sites.yale.edu/sites/default/files/BdM_A%20New%20Model%20for%20Predicting%20Policy%20ChoicesREvised.pdf
4. http://www.scirp.org/journal/PaperDownload.aspx?paperID=49058
5. The Predictability of Foreign Policies, The British EMU Policy, https://www.rug.nl/research/portal/files/3198774/13854.pdf
6. J. Velev, Python Code, https://github.com/jmckib/bdm-scholz-expected-utility-model.git

Answer 1

This is answered completely with a Python implementation by Jeremy Velev here: https://github.com/jmckib/bdm-scholz-expected-utility-model

The details are very hard to get right, look here for commentary and another Python implementation by David Masad here: http://davidmasad.com/blog/bdm-model-sensitivity-analysis/

These authors and Scholz et al note many issues with BDM including:

• BDM assumes that risk tolerances are the same for all agents
• Notations change constantly between publications and can be inconsistent, incomplete or entirely undefined
• No confidence interval is supplied with the prediction, and sensitivity analysis shows that the predictions for the examples given have a very wide confidence interval