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() / \

    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
            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)
            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)
            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
            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
                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]"

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
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.



Labor Pro EMU,100,75,40
Labor Eurosceptic,50,35,40
The Bank of England,10,50,60
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



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.

Output plot 1

Output plot 2

Output plot 3

Output plot 4

Output plot 5

  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
  • \$\begingroup\$ Folks I think he's asking about the results, as they compare to de Mesquita's published claims, not the whitespace in his code! David Masad has covered this topic here and also published Python code. \$\endgroup\$ Apr 20 '17 at 17:16
  • 1
    \$\begingroup\$ OP is trying to write a Python program to reproduce a claimed calculation result of Bueno De Mesquita (BDM). There is another, equally conscientious attempt to reproduce this calculation, in Python, by David Masad, "Replicating a replication of BDM". Masad also provides Python code, and also shows an approximately 20% divergence in the median score, starting from the same example and same inputs and same references. It is reasonable to conclude that BDM's paper does not provide adequate information to reproduce BDM's results. \$\endgroup\$ Apr 20 '17 at 23:22

First off, this:

Q = 1.0 ; T = 1.0

Should be expanded to this:

Q = 1.0
T = 1.0

And, if Q and T aren't constants, they should be renamed to q and t.

Secondly, you're missing whitespace in lots of places. For example, this:

def v(self,i,j,k):

Should be, again, expanded to this:

def v(self, i, j, k):

You should also have whitespace between mathematical operators as well, like this:

((1 + 2 * (2 - 6)) / 10) % 5

In addition to mathematical operators, you should also have whitespace between comparison operators as well, like this:

True != False and 50 <= (50 ** 2)

Adding whitespace in general where it's needed will help improve the overall readability of your code.

Your naming is, well, not the best. While this is scientific/mathematical code, I'd still recommend naming things better. For example, I have no idea what the function Ubi_i does, or what it's arguments, i, j, or ri do. Preferably, the names should be more descriptive of what purpose a variable/function/class serves.

A few bad names might include:

  • df
  • tmp1
  • tmp2
  • v
  • EU_i
  • l
  • res
  • eui
  • euj

On the topic of naming, Python has an official style guide, PEP8, and the style for naming is as follows.

  1. Variables/function arguments should be in snake_case, but if their value is constant, it should be in UPPER_SNAKE_CASE.
  2. Functions should be in snake_case.
  3. Classes should be in PascalCase.

When creating classes in any version of Python 2.x, you need to have the class explicitly inherit from object, like this:

class MyClass(object):

If you're using any version of Python 3.x or higher, this isn't required, and you can just do this instead:

class MyClass:

You have a couple of places where you're mixing single, '', and double quotes, "". Usually, it's best to be consistent.

You should also add some docstrings to your functions/classes. Docstrings should provide a description of what the function/class does, and it's arguments, if it has any. Here's an example of a docstring.

def my_func( ... ):
    Describe your function/class and it's
    arguments here.

Finally, in your run.py file, all the code, (except for the imports), should be wrapped in an if __name__ == "__main__": block to ensure that it properly runs. See this Stackoverflow question for more details.

  • 1
    \$\begingroup\$ On your first point, Q = T = 1.0 or Q, T = 1.0, 1.0 would also be acceptable, and a little neater. And your docstring doesn't quite meet PEP-0257. \$\endgroup\$
    – jonrsharpe
    Jul 18 '15 at 8:35
  1. Why did you write your own functions for mean and weighted_median? It looks like you are doing fairly standard things with them, so I would use implementations from numpy or scipy instead. (And calling the median of a vector a times w a "weighted median" of vector a makes some sense but it is really just the regular vanilla median of a*w.)

  2. Your Game class seems pretty monolithic. I'd either remove it entirely and just use global-level functions and scripts, or make more classes. For example it looks like Round could be a class, and maybe even Actor.

  3. I agree with all of Ethan Bierlein's comments too.

  4. You could make the binary logic statements in the big decison making loop a bit easier to read. For example, instead of duplicating np.abs(eui) < np.abs(euj) and np.abs(eui) > np.abs(euj) in every elif, you could pre-define a abs_eui_greater = np.abs(eui) > np.abs(euj) before the if and then just use that variable instead of repeating the binary logic. You could do something similar with the other repeated logic.

  5. In a few spots you are mixing numpy and native Python when you don't have to. For example, does tmp1 = np.sum(map(lambda (i,j,k): self.v(j,i,k), lpos)) really need np.sum() or would the python built-in sum() do the job? Or even better yet, if why not make the Salience, Ui_ij, etc. variables numpy arrays?


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.